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CHAPTER  1 
INTRODUCTION 

The  normal -mode  theory  gives  the  exact  solution  to  the  problem  of  sound 
propagation  in  a  range  -  independent  oceanic  wave -guide.  It  is  well  known, 
however,  that  the  ocean  has  some  range  dependence  due  to  changes  in  the  depth 
of  the  bottom  or  changes  in  the  acoustic  properties  of  the  environment  as 
shown  in  Figure  1.  From  all  the  acoustic  properties,  the  sound  speed  is  the 
most  range -variable  one  due  to  the  changes  in  temperature  and  pressure. 
Therefore,  a  reliable  range -dependent  method  must  be  created.  This 
range -dependent  method  must  allow  for  the  effects  of  bottom  interaction, 
including  shear  wave  from  the  solid  sediments  and  the  attenuation  coefficient 
of  the  compressional  and  shear  wave.  The  most  suitable  method  for  this 
purpose  is  the  coupled  normal -mode  method.  This  method  consists  in  dividing 
the  range -dependent  environment  into  range  - independent  segments,  solve  for 
the  eigenvalues  and  eigenfunctions  of  these  segments  using  the  normal-mode 
theory,  and  couple  these  modes  using  the  coupled  range  equation  and  the 
radial  boundary  conditions. 

The  coupled  normal -mode  method  was  implicitly  originated  by  Allan  Pierce 
in  the  late  60's  as  an  adiabatic -mode  theory,  with  eigenray  calculations  to 
estimate  the  coupling  coefficients,  to  simplify  the  solution  to  the 

1  -  4 

propagation  in  a  range -dependent  environment.  In  this  adiabatic  mode 

theory,  he  assumes  isovelocity  media  and  a  weak  coupling  between  the  natural 
modes  of  the  wedge -like  wave  guide.  He  found  that  compressional  waves 
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refract  into  the  basement  until  they  get  completely  attenuated,  and  his 
results  agreed  well  with  the  results  from  the  Parabolic  Equation  method. 
McDaniel  used  the  coupled  eigen-equations  to  calculate  the  energy  transferred 
between  normal  modes  as  a  result  from  bottom  scattering  of  the  ocean,  and  has 
shown  that  randomly  rough  sea  bed  layering  can  increase  the  transmission  loss 
depending  upon  the  degree  of  penetration  of  the  acoustic  field  into  the 

5  ~  8 

sediment.  Evans  modeled  the  axisymmetric  range -dependent  medium  as  N 

range- independent  segments  with  a  pressure-release  false  bottom  suitably  deep 

to  discretize  the  continuous  spectrum,  and  isodensity- isovelocity  layers  to 

simplify  the  solutions  of  the  linear  wave  equation  in  each  layer.  In 

this  stepwise  coupled-mode  method  he  solves  for  the  eigenvalues  and 

eigenfunctions  of  each  range- independent  segment  by  taking  into  account  the 

absorption  of  the  basement  to  avoid  the  reflected  energy  from  the  pressure - 

release  false  bottom.  Uberall's  group  has  used  a  similar  method  for  solving 

the  set  of  coupled  range  equations,  but  assuming  layers  of  linear  index  of 

1  2  w  2  C 

refraction  squared.  This  gives  a  better  approximation  of  the  sound 

speed  profile  and  the  eigenfunctions  are  Airy  functions.  However,  the 
effects  of  compressional  attenuation  of  all  water  and  sediment  layers  is  not 
included  in  any  of  the  existing  models. 

Many  theoretical21  37  and  experimental3  8  5  4  scientists  have  shown  the 
importance  of  attenuation  and  shear  waves  in  the  calculation  of  the 
transmission  loss  and  the  study  of  sound  propagation  in  the  ocean.  Low 
frequency  compressional  and  shear  waves  are  less  attenuated  in  all  layers,  so 
they  can  penetrate  large  depths  and  distances  of  the  ocean.  In  view  of  the 
current  interest  in  low  frequency,  active  and  passive  sonar  systems,  it  is 
desirable  to  extend  the  Normal  Mode  theory  (which  is  an  exact  method 
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throughout,  but  is  most  easily  applicable  to  the  low  frequency  case  where  the 
number  of  modes  is  small)  to  the  situation  of  an  environment  containing  range 
dependence,  absorption  in  all  the  layers  of  the  media,  and  birefringence 
effects  due  to  a  shear  speed  differing  from  compressional  speed  on  the  solid 
bottom  sediments.  Assuming  high  frequencies,  the  very  graphic  techniques  of 
acoustic  ray  tracing  have  become  the  standard  tools  for  transmission  studies 
in  the  deep  ocean.  However,  the  more  exact,  but  less  graphic,  normal  mode 
theory  is  particularly  suited  to  low  frequency  propagation  in  a  shallow 
ocean.  Normal  mode  calculations  at  higher  frequencies  are  also  possible 
whenever  computer  time  and  memory  are  available. 

The  formalism  to  be  developed  will  be  based  on  an  existing  normal  mode 
propagation  model  for  a  range  dependent  environment.12  20  This  model  admits 
general  depth  and  range  dependent  sound  velocity  profiles  and  an  arbitrary 
but  gradual  range  dependence  of  a  layered  ocean  bottom  containing  sound 
velocity  gradients.  Sound  penetration  into  the  bottom  is  accounted  for  by 
considering  not  only  the  modes  trapped  in  the  water  column  but  also  a 
sufficiently  large  number  of  radiating  modes.  The  set  of  radiating  modes, 
which  is  in  principle  continuous,  is  discretized  by  a  basement  layer  of  large 
but  finite  thickness.  For  a  range -dependent  environment,  a  set  of  coupled 
range  equations  (one  for  each  mode)  is  obtained  which  is  solved  by 
diagonalization  procedures.  The  algorithm  includes  absorption  in  the 
basement  as  a  first  order  perturbation  to  the  problem  and  it  ignores  the 
shear  waves  created  in  the  solid  sediments  of  the  bottom. 

It  is  intended  here  to  upgrade  this  existing  model  first  by  taking  into 
account  the  compressional  attenuation  coefficient  of  all  the  layers  of  the 
medium  and  to  provide  a  detailed  derivation  of  the  coupled-mode  theory. 
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Eventually,  in  a  next  report,  the  generation  and  propagation  of  shear  waves, 
with  the  effects  of  its  absorption,  in  the  bottom  sediments  will  be  included 
to  perform  quantitative  studies  on  the  effects  of  absorption  and  shear  waves 
in  the  propagation  of  sound  in  the  water  and  bottom  of  a  range  -  dependent 
stratified  ocean.  It  has  been  shown  that  the  exact  treatment  of  shear  waves 
in  models  such  as  the  parabolic  equation  (PE)  approximation  is  not 

....  55.56 

possible . 

Absorption  effects  will  be  accounted  for  as  a  first  and  second  order 
perturbation  to  the  simpler  problem  of  propagation  in  a  non-absorbing 
environment  instead  of  using  complex  sound  speed  profiles,  leading  to  complex 
eigenvalues  and  eigenfunctions  since  this  will  necessitate  a  transition  from 
real  to  complex  Airy  functions  as  the  depth  function  solutions  in  each 
layer.57  Including  absorptive  effects  will  not  only  make  it  possible  to 
simulate  more  realistic  ocean  environments  but  it  will  also  lead  to  a  more 
satisfactory  procedure  for  discretizing  the  continuous  distribution  of 
radiating  modes.58 
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CHAPTER  2 
THEORY 

2.1  DERIVATION  OF  THE  WAVE  EQUATION  FOR  FLUID  LAYERS 

The  wave  equation  is  a  mathematical  description  of  the  reaction  of  the 
media  due  to  a  disturbance  from  an  external  force  caused  by  a  source  or 

sources.  The  media  can  well  be  a  gas,  liquid,  or  a  solid,  and  the  source  may 

be  electromagnetic  or  mechanical.  In  this  report,  the  mechanical  (acoustic) 
propagation  of  the  disturbance  in  a  fluid  media  is  treated.  The  fluid  media 
is  the  ocean  environment  modeled  as  a  horizontally  stratified  acoustic 
waveguide  were  the  surface  is  treated  as  a  pressure-release  (resilient) 
boundary  and  the  bottom  is  taken  as  sediment  layers  with  variable  sound 
speed,  density,  and  attenuation  coefficient.  The  effects  of  bottom  shear 
waves  are  neglected  even  though  their  importance  is  well  understood  and 

documented  in  books  as  the  ones  in  References  59  to  62.  Shear  wave  and  its 

attenuation  coefficient  will  be  included  to  these  calculations  in  a  future 
report . 

The  disturbance  created  by  the  acoustic  source  may  be  expressed  as  a 
change  in  the  total  pressure,  relative  to  the  undisturbed  pressure,  as  a 
function  of  the  density  fluctuation  created  by  this  external  force.  If  the 
density  fluctuation  is  much  smaller  than  the  undisturbed  density  of  the 
environment,  then  the  total  disturbed  pressure  may  be  expanded  in  the 
following  Taylor  series: 
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P(p)  “  PQ+ 


3P 


ap 


<?•/.„)  +  T 


f  3  P  1 


Bp 


(P-P0)‘  +  ... 


(1) 


o  r  "  o 

where  the  partial  derivatives  are  constants  determined  for  the  adiabatic 

compression  and  expansion  of  the  fluid  about  its  equilibrium  density  p  ,  the 
equilibrium  pressure  is  P  ,  and  the  instantaneous  total  density  is  p. 


If  the  magnitude  of  the  condensation  is  much  smaller  than  unity,  i.e. 


s  ■  (p  -  PQ)/P0  ■  P_/PQ  (2) 

then  then  first  two  terms  in  the  Taylor  expansion  are  of  greatest 
contribution  and  an  acoustic  pressure  caused  by  the  disturbance  may  be 
defined  as 


p  -  P(p)  -  Pq  a 


<9P 


l  dP  J 


(3) 


where  by  thermodynamic  arguments  it  is  found  that  in  an  adiabatic  media  the 
sound  speed  is  given  by 


and  the  adiabatic  bulk  modulus  is  given  by 

B  -  p  c2 
o 

therefore,  the  acoustic  pressure  is  simplified  to 

2 

p  =>  c  p 


(*> 

(5) 

(6) 


which  is  called  the  equation  of  state. 

An  equation  for  the  motion  of  the  particles  in  the  fluid  is  also 
necessary  for  the  proper  environmental  description.  Consider  an 
infinitesimal  cubic  volume  in  the  media  where  the  disturbance  is  taking  place 
as  shown  in  Figure  2(a)  for  the  one -dimensional  derivation  in  cartesian 
coordinates.  Equating  forces  in  a  continuous  medium  gives 
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S' 


F  +  A  [P(x)  -  P(x+dx) 1  -  J-(p  A  dx  jg)  ( 

•xtemal  at  Ot 

where  the  external  force  is  the  disturbance  created  by  the  sound  source  and 

can  be  written  in  terms  of  a  "force  density"  with  the  expression 


F  -  X  A  dx 

•xtemal  • 


and  using  the  definition  of  a  derivative 

8?  P(x+dx) -P(x) 


dx 


dx 


Equation  (7)  becomes 


x.  -  I  “ v,  k(V>  +  St<V> 


which  in  three  dimensions  is  given  by 


Xe-  VP  -  V  (V-pV)  +  |^(pV) 


(8) 

(9) 

(10) 


(11) 


where  the  density  is  a  time  and  space  function  and  the  equation  is  in  a 

non-linear  form.  The  total  pressure  is  P  and  the  total  instantaneous 

— ♦ 

particle  velocity  is  V.  Dividing  the  instantaneous  density,  pressure,  and 
particle  velocity  into  an  undisturbed  part  and  an  acoustic  part  Equation  (11) 
simplifies  to  the  linear  form 

3p  +  -  xt  (12) 

where  the  undisturbed  density  of  the  medium  is  p  ,  the  acoustic  pressure  is 

—¥ 

p,  and  the  particle  velocity  is  v. 

Since  the  fluid  of  interest  is  continuous  throughout  the  infinitesimal 

volume,  Figure  2(b)  will  be  helpful  in  deriving  a  continuity  equation  under 

the  argument  that  the  continuous  mass  going  into  the  volume,  p(x)  A  V  (x)  dt, 

must  be  the  same  as  the  quantity  going  out,  p(x+dx)  A  V  (x+dx)  dt .  There  may 

be  an  change  in  mass  inside  the  volume  due  to  the  compressibility  of  the 

fluid,  ^  A  dx  dt,  and  there  may  exist  a  source  of  mass  inside  the  volume 
dt 

represented  by  Q  A  dx  dt.  Taking  the  definition  of  the  derivative  in 
Equation  (9)  gives 
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*  4-(pV)  A  dx  dt  -  ^  A  dx  dt  +  Q  A  dx  dt 
dx  3 1 


which  is  rewritten  in  three  dimensions  as 


V  •  (pV)  +  ||  +  Q  -  0 


(13) 


(14) 


or  in  a  linearized  form  as 

■  5  +  ff  - 0  <15> 

63 

where  Q-0  for  the  sources  of  interest. 

Substituting  Equation  (6)  into  Equation  (15)  for  the  acoustic  density, 
taking  its  partial  derivative  in  time,  and  dropping  the  subscript  0  of  the 
undisturbed  density  of  the  medium,  gives 


d  -*  -*  -id 
P  fr(V-v)  +  c  p  -  0 

at  St2 

and  taking  the  divergence  of  Equation  (12)  yields 

pV(l/p  Vp)  +  p  |^(V-v)  -  p  V-(X/p) 


(16) 


(17) 


which  subtracted  from  Equation  (16)  provides  the  inhomogeneous  wave  equation 

pV-(l/p  Vp)  +  k2p  -  p  V • ( 1/p  VU)  (18) 

where  the  external  force  has  been  written  in  terms  of  an  external  potential 
energy,  time  harmonic  behavior  has  been  assumed  where  k-«/c ,  and  the 
undisturbed  density  of  the  fluid  is  taken  as  space  dependent.  This  equation 
can  also  be  written  as 


v2p  +  k2p  -  p_1(Vp).(Vp)  -  V2u  -  p'J(Vp)  •  (Vt/) 


which  is  simplified  under  the  change  of  variables 

p  -  /p  n 

and 

U  -  /p  v 

to  obtain 


6* 


V2n  +  (k2+AT2)n  -  v2t/  +  K2 V 


(19) 


(20) 


(21) 


(22) 


where 


8 


NSWC  TR  88-166 


7i  ■  ?C  7  fy2-  <23> 

If  the  density  is  taken  as  a  linear  function  of  depth,  then  Equation  (23) 
simplifies.  However,  the  inhomogeneous  equation  to  solve  also  has  a  depth 
dependent  wavenumber  to  worry  about  due  to  the  depth  dependence  of  the  sound 
speed.  The  changes  in  density  with  depth  hardly  occurs  compared  to  the 
changes  in  sound  speed.  It  is  concluded  that,  for  simplicity,  the  ocean 
environment  can  be  divided  into  horizontal  layers  with  constant  density.  It 
is  understood  that  the  bottom  sediments  may  have  layers  of  large  density 
gradients.  In  this  case,  Equation  (22)  must  be  solved.  In  this  report, 
however,  layers  of  constant  density  are  assumed.  Then  Equation  (18)  becomes 

V2p  +  k2p  -  V2U  (24) 

which,  in  an  unbounded  medium,  has  a  general  homogeneous  solution  consisting 
of  an  outgoing  and  an  incoming  wave,  and  an  inhomogeneous  solution  caused  by 
the  external  force.  Since  a  sound  source  in  a  fluid  can  only  produce  a 
scalar  potential  (no  shear  waves),  the  curl  of  Equation  (12)  takes  us  to  the 
property  that 


— ►  — *  _ 

Vxv  -  constant  ■  vorticity  (25) 

the  vorticity  in  the  media  does  not  change.  Therefore,  if  initially  their 

has  been  no  rotational  component  of  the  particle  velocity  then  the  vorticity 

will  always  be  null  and  this  particle  velocity  can  be  written  in  terms  of  a 

velocity  potential 


v  ■  Vtp 

which  substituted  back  in  Equation  (12)  gives 


and  assuming  harmonic  time  dependence  yields 

p  -  -t up  (f 


(26) 


(27) 


(28) 
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or 


-tVp 


up 


(29) 


Substitution  of  Equation  (28)  into  Equation  (24)  provides  the  inhomogeneous 
Helmholtz  equation 


V2cp  +  k  2<p - —  V2U  (30) 

u  p 

which  must  be  solved  for  the  velocity  potential.  If  the  medium  is  bounded, 
then  the  solution  must  satisfy  the  appropriate  boundary  conditions. 

The  conservation  of  energy  is  obtained  by  dotting  Equation  (12)  with  the 
particle  velocity  and  substituting  the  continuity  equation,  Equation  (15) , 
giving  the  conservation  law 

Yz  +  V-I  -  0  (31) 


where 

J?  -  \  P  v2  +  (32) 

Z  2  pc 

is  the  acoustic  energy  density,  and 

I  -  p  v  (33) 

is  the  acoustic  energy  flux  or  acoustic  intensity.  Integrating  Equation  (31) 
throughout  a  volume  in  the  fluid  medium  provides  the  power 

n-^Jl?dV--^I.ndS  (34) 

V  s 

in  terms  of  a  closed  surface  integral  around  the  volume  where  all  the  energy 
is  contained. 

To  obtain  these  important  measurable  quantities,  it  is  necessary  to 
solve  Equation  (30)  for  the  velocity  potential.  To  solve  this  equation  by 
separation  of  variables,  it  has  been  proven  that  the  sound  speed  must  be 
function  of  one  variable.  This  variable  is  taken  to  be  the  vertical 
direction  since  temperature  and  the  total  pressure  of  the  ocean  highly 
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dependence  on  depth.  A  simplified  version  of  Wilson's  formula  for  the  sound 

speed  as  a  function  of  temperature,  salinity,  and  depth  is  given  by 

c(z)  -  1492.9  +  3(T- 10)  -  6xlO*3(T- 10)2  -  4xl0'Z(T- 18)2  +  1.2(S-35)  - 

10*2(T-18) ( S -  3  5 )  +  z/61  (35) 

where  the  temperature  T  is  in  Celsius,  the  salinity  S  is  in  parts  per 
thousand,  and  the  depth  z  is  in  meters.  The  formula  is  accurate  to  0.1  m/s 
for  a  temperature  less  than  20°C  and  for  depths  less  than  8.0  kilometers.65 
Range  dependence  of  the  sound  speed,  density,  and  position  of  the  boundaries 
are  taken  care  of  by  the  coupled  normal-mode  method  in  Sections  2.6,  2.7,  and 
2.8  of  this  report. 
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2.2  SOLVING  THE  WAVE  EQUATION 


It  has  been  found  that,  in  the  liquid  layers,  the  inhomogeneous 


Helmholtz  equation 


[V2  +  k2(r) }  <p - —  V2U 

Ct)  P 


must  be  solved,  where  the  sound  speed  is  assumed  to  be  a  function  of  depth 


only . 


In  the  case  of  a  point  source,  Equation  (36)  may  be  written  as 


[V2  +  k2<0]  *  -  -  S  6(r  -  r  ) 

O 


where  S  contains  all  the  constants  in  the  inhomogeneous  term  and  is  usually 
called  the  source  strength.  For  simplicity  and  without  loss  of 
generalization  we  may  set  S  -  1  corresponding  to  a  unitary  source  strength, 
then  in  cylindrical  coordinates  the  equation  becomes 


1  d  f  d  1  ,  d2  ,2,  .1  ,  .  1  . ,  .  , 

;  Si  r  3i  +7T  +  k<z)  -  ■  Ri  5<r)  ‘ 


(z  -  z  ) 

O 


where  the  source  is  at  r  -  0,  z  -  z  ,  and  the  wavenumber  is  taken  to  be  depth 

5  8 

dependent  only  since  the  sound  speed  must  be  a  function  of  one  variable. 
Substituting  the  Fourier-Bessel  transformation, 


the  closure  relation, 


the  Bessel  equation, 


00 

<p(r,z)  -  f  u(k,z)  J  (kr)  k  dk, 
oJ 

u(k,z)  -  |  ip( r,z)  J  (kr)  r  dr, 
oJ 

00 

6 (r  -  r')  -  r  f  J  (kr)  J  (kr')  k  dk, 
Jo  0 
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k2r2  JJJ(kr)  +  kr  J'(kr)  +  k2r2  JQ(kr)  -  0,  (42) 

and 

o  r  ~  > 

%-  r  %-  J  (kr)  -  k  J'(kr)  +  k2r  J’’(kr)  -  -k2r  J  (kr)  (43) 

3r[  5r  j  o  o  o  o 

into  Equation  (38)  leads  us  from  a  partial  differential  equation  to  the 

ordinary  differential  equation 

-  ,2  -1  S(z  -  Z  ) 

— -  +  k2(z)  -  k2  u(k,z)  - - 2-  °  (44) 

-  dz  -< 

2 

where  k  is  the  eigenvalue  and  u(k,z)  is  the  eigenfunction  of  the 
inhomogeneous  equation.  After  solving  for  u(k,z),  we  can  transform  back 
using  Equation  (39)  to  obtain  ip(r,z).  The  solution  of  this  inhomogeneous 
equation  can  be  written  as  the  sum  of  the  homogeneous  solution  and  the 
particular  or  transient  solution.  The  generalized  homogeneous  solution  can 
be  used  as  the  solution  of  Equation  (38).  The  solution  of  Equation  (39), 
however,  is  different  since  the  shear  speed  and  compressional  speed  of  a 
sediment  are  usually  unequal . 

Since  the  bottom  of  the  oceanic  waveguide  is  not  rigid  nor  resilient, 
its  energy  spectrum  will  contain  discrete  modes  and  a  continuous  distribution 
which  represent  the  energy  that  radiates  to  infinity.  The  constant 
wavenumber  of  the  bottom  basement  represents  the  dividing  point  between  the 

2 

continuous  and  discrete  modes.  It  is  possible,  however,  to  discretize  the  k 

spectrum  by  adding  to  the  problem  a  false  resilient  boundary.  In  this  case, 

the  homogeneous  equation  to  solve  becomes 

-  ,2 

—  +  k2(z)  -  k2  u  (z)  -  0  (45) 

,  2  n  n 

L  dz  J 

where  n  -  1,2,..., N  is  the  mode  index.  This  false  bottom  must  be  deep  enough 
to  make  the  discretized  radiative  modes  very  close  to  each  other  and  better 
simulate  the  continuous  spectrum.  Also  it  has  to  be  highly  absorptive  to 
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avoid  reflections  of  the  modes  with  the  false  bottom.  Since  the  Earth  is  not 
flat,  a  false  resilient  bottom  satisfying  the  above  conditions  is  an 
approximation  as  good  as  the  assumption  of  a  semi -  inf inite  basement  or  that 
of  a  rigid  false  bottom. 

From  the  continuity  of  pressure  in  the  liquid  layers,  from  Section  2.3, 
we  shall  set  the  function  /p(z)  u^(z)  as  the  orthonormal  eigenfunctions, 

|  p(z)  u*(z)  u  (z)  dz  -  S  (46) 

J  n  ro  nm 

with  the  closure  relation 


6(z  -  zq)  -  p 


and  the  eigenfunction  is  given  by 


n 

(Z  )  )  U  (z)  U*(z  )  , 
0  _ (  n  n  0 


(47) 


n  *  1 


u(k,z)  -  )  a  (k)  u  (z)  . 

L _ .  n  n 


(48) 


n  -  1 


The  inhomogeneous  term  is  taken  into  account  if  we  substitute  the  homogeneous 
solution  in  the  inhomogeneous  equation.  The  substitution  gives 

r  d2 


+  k2(z)  -  k2 


dz4 


an(k) u, 


.  .  '<zo> 

n(Z) - 2tT 


)  u  (z)  u 

L _ .  n  n 


(Z0)  (49) 


n  -  1 


n  »  1 


and  substituting  the  homogeneous  equation  we  get 

N 


Eh 


p(z  ) 

(k)  (k2  -  k2)  -  -5-2-  u‘(z  ) 

n  Z7T  n  0 


U  (z)  -  0 


(50) 


n  -  1 


and  to  satisfy  the  equation  we  must  set  the  terms  inside  the  brackets  to  zero 
leading  to 


p(z  )  u  (z  ) 

/i  \  0  n  0 

a  (k) - -z -  - 

n  /7T 


k2  -  k2 


(51) 


which  substituted  into  Equation  (52)  gives 
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u(k,z)  - 


2tt 


r — .  u  (z  )  u  (z) 

\  n  0  n 

l—.  u2_  1* 


(52) 


k  -  k 


n  *  1  n 

and  this  substituted  into  Equation  (44)  gives  the  scalar  potential 

N 


p(zo)  V-’  * 

cp(r,z) - tt—  )  u  (z  )  u  (z) 

ZJT  [_ _ _  n  0  n 


f®  j^(kr)  k  dk 


n  “  1 


OJ 


k2  -  k2 

n 


(53) 


where  the  integral  of  this  equation  is  better  solved  by  contour  integration 
and  we  may  define  it  as 


Vr)  -  2 


H(1>(kr)  +  H{2>(kr) 
o  o 


0j 


,2  .2 

k  -  k 


k  dk. 


(54) 


A  property  of  the  eigenvalues  of  the  problem  is  that  these  have  a  very 
small  imaginary  part  compared  to  the  real  part.  Also,  both  parts  of  the 
eigenvalues  are  positive.  This  is  because  we  consider  only  outgoing  waves 
from  the  source.  To  solve  the  integral  we  can  consider  contour  integration 
in  the  first  quadrant  of  the  complex  k-plane.  Consider  the  contour  path 
displayed  in  Figure  3,  where 


I  (r)  -  ^ 

ni 


H(1,(kr) 
o 

2  2 
k  -  k 

C  n 

1  i 


H<2>(kr) 

o 


k  dk,  i  -  1,2,3 


(55) 


i  2  i  2 

k  -  k 

C  n 

2  i 


k  dk,  i  -  4,5,6 


and,  by  this  definition,  the  integral  we  want  to  solve  is, 

l 


I  (r) 


I  (r)  +  I  (r) 

nl  n4 


6  3 


By  Jordan's  lemma  we  have  that 


I  (r)  -  I  (r)  -  0 

n2  n5 

also  the  integrals  in  Equation  (55)  show  that 


(56) 


(57) 
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I  (r)  +  I  (r)  -  0 

n3  n6 


which  means  that  we  may  write 


r  3  6  -|  6 

I  (r)  -  |  V  I  (r)  +  V  I  (r)  T  I  .(r)  (5? 

n  2  U  ni  L.  ni  2  L.  ni 

L  i  - 1  i - *  J  i-1 

where  only  I  (r)  and  I  (r)  contribute  to  the  sum. 

nl  n* 

Given  that  the  singularities,  k  -  k  ,  are  located  in  the  upper  contour 


we  get  that 


6 

y  I  (r)  -  O 

l—  ni 

i«4  J 


„<  2  )  . 
H  (kr) 
o 

.2  ,2 
k  -  k 


k  dk  -  0 


l  Vr>  - 1 

i-l  J 


*  H( 1 5 (kr) 

1  O  — - -  k  dk  -  rc  L  H(1)(k  r) 

,2,2  O  n 

J  k  -  k 

C  n 

1 


by  calculus  of  residues.  Substituting  Equations  (60)  and  (61)  in  Equation 
(59)  and  this  one  into  Equation  (53)  gives 


<p(r,z)  -  ^  p( z  )  V  u*(zn)  u  (z)  Hln(k  r) 
**  0  Z _ _  n  0  n  On 


where  the  eigenfunctions  u  (z),  and  eigenvalues  k  ,  satisfy  Equation  (45). 

n  n 

Note  that  the  solution  can  be  written  as  a  separation  of  variables.  Figure  4 
shows  a  qualitative  picture  of  the  shape  of  the  modes  for  several  sound  speed 
profiles . 

Now  we  must  solve  the  characteristic  equation  for  the  compressional 
waves,  Equation  (45).  We  may  start  by  dividing  the  sound  speed  profile  into 
layers  where  the  squared  of  the  index  of  refraction  is  a  linear  function  of 
depth,  or 

n2(z)-az+b  (63) 
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where  k(z)  -  u  n(z)  and  Figure  5  gives  the  geometry  to  be  used  for  this 
investigation.  To  determine  a  and  b  we  let  the  sound  speed  at  the  top  of  the 
layer  to  be  c  and  that  of  the  bottom  to  be  c  .  Substituting  into  our  linear 
equation  gives 

“  a  zt  +  b,  (64) 

c 

t 

and 


—  -  a  z  +  b , 

2  b 


(65) 


which  solved  for  a  and  b  results  in 


2  2 

c  -  c 
t _ b_ 

.22 
€  C  C 
t  b 


(66) 


and 


b  - 


,  2  2\ 

z  (c  -  c  ) 

t  t  b 

,22 
l  C  C 
t  b 


(67) 


where  l  -  z  -  z  is  the  thickness  of  the  layer.  If  we  define  the  indices  of 

b  t 

2  2  2  2  2  2 

refraction  n  -  1/c  ,  n  -  1/c  ,  and  a  -  (n  -  n  )/l  ,  then  we  have 

t  t  b  b  t.  b 


,2,  .  2 
k  (z)  -  u 


n  +  a  (z  -  z) 

t  t 


k‘  +  S  (z  -  z  ) 


where  S  -  Ak  /Az ,  and  which  substituted  into  our  eigenequation  gives 

j2 


dz 


-  u  (z)  + 

2  n 


k2  +  S  (z  -  z  )  -  k2 
t  t  n 


u  (z)  -  0 . 


Define 


f  (z) 


-  S 


-2/3 


k2  +  S  (z  -  z  )  -  k2 

t  t  n 


and  square  its  derivative  to  obtain 

- - |  s  j 2/3  — 

2  1  1  2 
dz  d  r 

which  substituted  into  the  new  eigenequation  gives 


(68) 


(69) 


(70) 


(71) 
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-  f  u  (O  -  0.  (72) 

L  df2  -i  n 

The  solutions  of  this  differential  equation  are  the  Bessel  functions  of  order 
1/3,  or  more  commonly  known  as  the  Airy  functions,  i.e. 

u  (f)  -  a  Ai(?)  +  b  Bi(0  .  (73) 

n  n  n 

The  behavior  of  these  functions  and  their  derivative  is  shown  in  Figure  4. 

Now  that  the  general  solutions  are  found,  we  must  match  the  solutions  at  each 
boundary  with  the  appropriate  boundary  conditions  in  order  to  find  the 
unknown  coefficients  and  eigenvalues  for  each  mode. 
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2.3  DEPTH  BOUNDARY  CONDITIONS 

2.3.1  Boundary  Between  Fluid  Layers 

The  boundary  conditions  for  the  interface  between  two  fluid  layers  are 
obtained  when  an  infinitesimal  cylindrical  volume  is  modeled  across  this 
boundary.  In  this  case  we  have  to  satisfy  two  boundary  conditions. 


2. 3. 1.1  Continuity  of  the  Normal  Particle  Velocity.  The  volume  integration 
of  Equation  (15)  in  this  infinitesimal  cylinder  provides  the  expression 

p£v-ndA--|^JpdAAx  (74) 

where  making  Ax  — )  0,  the  right  hand  side  of  the  equation  vanishes  and  the 
surface  integral  yields  the  boundary  condition 

A  A 

-*■  -f 

(75) 


v  •  n  -  v  •  n 

2  s  1 


This  boundary  condition  is  expressed  as  v  •  n  -  continuous,  and  in  the 
case  of  stratified  layers  we  may  write  n  -  z  to  obtain  the  boundary 


condition, 


or  using  Equation  (62)  we  get 


dip 

v  -  —  -  continuous 
z  dz 


du 

n 

d T 


-  continuous. 


(76) 


(77) 


2. 3. 1.2  Continuity  of  the  Pressure.  Assuming  that  there  is  no  source  in 
the  infinitesimal  volume  of  this  cylinder,  the  volume  integration  of  Equation 
(12)  gives 

-  (£  p  n  dA  -  p  [vdAAx  (78) 
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and  letting  Ax  — >  0  yields  the  pressure  boundary  condition 

piL“  p2L 

where  substituting  Equation  (62)  into  Equation  (28)  and  this  one  into 
Equation  (79)  gives 

p  u  —  continuous. 


(79) 


(80) 


2.3.2  The  Resilient  Boundary  of  a  Fluid 

As  a  very  good  approximation,  we  may  consider  this  boundary  as  pressure 
release  for  acoustic  waves  in  the  liquid  layer.  Therefore,  the  only  boundary 
condition  is  that  the  pressure  vanishes  at  this  boundary,  i.e. 

u  |  -  0.  (81) 

n  z 

0 


2.3.3  Down-Layer  Matching  Algorithm 

Now  that  the  boundary  conditions  are  determined,  a  matching  algorithm  is 
next  to  by  determined.  The  normalization  constant  for  the  eigenfunctions  are 
given  by  the  expression 

J  ♦  1  Z 
r - -  -j  +  1  A 

N~2-  )  p  j  |u  .  |Zdz  (82) 

n  / .  J  J  nj  1 

-  z 

j  m  1  j 

A 

where  u  is  the  unnormalized  eigenfunction  and 
nj 


u  (z)«  N  u  (z) 

nj  n  nj 

is  the  normalized  one.  For  simplicity  we  will  make  some  further 

Ai.  *  Ai [ C  ,  (z  )  ] 

U  m  j 


a  * 

j 


dz 


Ai'  « 

i  j 


d  Ai 
dC  . 


? 

m 


(2  ) 

J 

-  sgn (5  ) 
j 


s 

J 


1/3 


(83) 

definitions : 

(84a) 

(84b) 

(84c) 
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and 


P  -  k2(z  )  (84d) 

j  j  j 


*>  -  (P  -  k2)/S  (84e) 

nj  j  n  j 

which  converts  Equation  (70)  to  the  form 

f  (z )-a  (z+b  ).  (85) 

nj  j  nj 

The  layering  subscript  runs  from  j-1  at  the  resilient  surface,  to  j-J+l»B  at 
the  upper  boundary  of  the  basement  as  shown  in  Figure  5. 

From  the  resilient  boundary  condition  at  the  surface  of  the  ocean, 
Equation  (81),  we  get 


u  (z  )  -  a  Ai  +  b  Bi  -  0 

nl  1  nl  11  nl  11 


(86) 


and  if  we  write 


a  *  T  Bi  , 

nl  D1  11 


b  -  -  T  Ai 

nl  D1  11 


then  the  solution  to  the  first  layer  is 


(87) 


u  (z)  -  T  [Bi  Ai(f  )  -  Ai  Bi(f  )]  (88) 

nl  D1 L  11  nl  11  nl  J  v 

where  we  can  arbitrarily  set  T  ■  sgn(Si)  to  ensure  a  positive  slope  of  the 
eigenfunction  at  the  top  layer  of  the  media. 

The  value  of  the  derivative  of  this  function  at  the  surface  simplifies 


to 


u*(zf) 

n  1  1 


T  [Bi  Ai’  -  Ai  Bi'  ; 

D1  11  11  11  11J 


where  the  term  inside  the  brackets  is  the  Wronskian  relation 

W  -  Ai(C)  Bi'(f)  -  Bi(C)  Ai' (D  -  1/n 


56 


(89) 


(90) 


then 


u' 

n  1 


(zi) 


T  ?r 


D1 


(91) 
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Going  back  to  Equation  (73)  for  the  boundary  at  z we  get 


u  (z  )  -  a  Ai  +  b  Bi 

nl  2  nl  12  nl  12 


u  (z  )  -  a  Ai  +  b  Bi 

n  2  2  n  2  22  n2  22 


du 

i 

dz 

A 

du 

i 

dz 


n  1 


n  2 


-  a  (a  Ai  +  b  Bi'  ) 
Z-Z  1  nl  12  nl  12 
2 


-  a  (a  Ai'  +  b  Bi'  ) 
Z-Z  2  n  2  22  n2  22 

2 


and  using  the  liquid  boundary  conditions  we  obtain 


and 


or 


p  (a  Ai  +  b  Bi  )  -  p  (a  Ai  +  b  Bi  ) 

1  nl  12  nl  12  2  n2  22  n2  22 


a  (a  Ai'  +  b  Bi'  )  -  a  (a  Ai'  +  b  Bi'  ) 

1  nl  12  nl  12  2  n2  22  n2  22 


R  u  (z  )  -  a  Ai  +  b  Bi 
D  2  n  1  2  n  2  22  n2  22 


(92a) 

(92b) 

(92c) 

(92d) 

(93a) 

(93b) 

(94a)  , 


and 


where  R 


Dj 


T  u’  (z  )  -  a  Ai'  +  b  Bi ' 

D  2  nl  2  n  2  22  n2  22 

p  /p  ,  and  T  -  a  /a  ,  where  j-2,3 . J. 

J-l  J  Dj  j-1  j 


In  matrix  form  we  write, 


R  u  (z  ) 

D2  nl  2 

A 

T  u'  (z  ) 

02  nl  2 


and  solving  for  a  and  b  we  have 

°  n2  n2 


Ai 


22 


Bi 


22 


Ai'  Bi' 
22  22 


1 

“ 

a 

n2 

b 

n2 

- 

a 

n2 

—  n 

b 

n2j 

R  Bi'  u  (z  ) 

D2  2  2  nl  2 


T  Bi  u'  (z  ) 

02  22  nl  2 


T  Ai  u'  (z  )  -  R  Ai'  u  (z  ) 

02  22  nl  2  D2  22  n  1  2 


or  in  general 


(94b) 


(95) 


(96) 


L 

— j 

L 

j 

r 

_ 

A 

A 

a 

nj 

-  n 

R 

Dj 

Bi' 

j  J 

u  .  (z  ) 
nj-  1  J 

A 

-  T 

D  j 

Bi 

j  j 

u' 

nj- 

A 

i(zj) 

(97) 

b 

nj 

T 

Dj 

Ai 

j  j 

U' 

nj-  1  j 

-  R 

D  j 

Ai' 
j  j 

u 

nj- 

1(V 

where  the  values  of  the  eigenfunction  and  its  derivative  at  the  bottom  of  any 
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layer  have  Co  be  determined  before  the  coefficients  of  the  next  layer  are 

sought.  This  solution  is  good  for  "normal -range"  values  of  f  . 

In  the  case  of  large  positive  values  of  f  we  can  substitute  a  and 

nj  nj 

b  into  the  general  solution  and  its  derivative  to  obtain  the  matrix 
nj 

equation 


u  (z) 
nj 

A 

U'.(Z) 

nj 


C  C  R  u  (z  ) 

11  12  Dj  nj-  1  j 

A 

C  C  T  u'  (z  ) 

21  22  Dj  nj-1  j 


where 


C  (T  )  -  Bi' .  Ai(f  )  -  Ai'  Bi(f  ) 

11  n  j  J  J  n  j  j  J  n  j 

c,,cr  .)  -  Ai  .  Bi(f  .)  -  Bi.  Ai(r  ) 

12  n  j  jj  nj  jj  nj 

C  (f  .)  -  Bi-  .  Ai'  (ST  .)  -  Ai ;  Bi '  ( f  .) 

2  1  n  j  jj  nj  jj  nj 

c,  (f  )  -  Ai  .  Bi'(f  .)  -  Bi. . Ai ' (f  .) 

22nj  jj  nj  jj  nj 

and  the  Airy  functions  nearly  cancel  each  other  for  large  f 


(99a) 

(99b) 

(99c) 

(99d) 


Now  that  a  recurrence  relationship  has  been  found  for  calculating  the 
eigenfunction  at  the  layer  j  after  calculating  it  at  the  layer  j - 1 ,  we  are 
left  with  the  boundary  conditions  at  the  basement.  This  basement  will  be 
taken  as  a  layer  with  constant  acoustic  parameters  since  we  will  make  it 
deeper  than  the  region  of  interest.  This  layer  is  created  for  the  only 
purpose  of  discretizing  the  continuous  wavenumber  spectrum.  However,  the 
eigenfunction  at  this  layer  is  exponential  in  shape  when  the  mode  is  trapped 
and  oscillatory  when  the  mode  is  a  radiating  one,  hence  we  must  divide  the 
solutions  into  two  cases. 


2. 3. 3.1  Trapped  modes  (k>k_).  The  eigenfunction  at  the  basement  is  given 

- - -  n  B 


a  -7  z 

u  (z)  -  a  e  +  b  e 

n  B  n  B  n  B 


(100) 
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2  2  2 

where  7  -k  -k^.  We  use  the  boundary  condition  at  z-z^  (resilient)  to  get 


b  -  -  a  e 

n  B  n  B 


-27  2 
n  F 


which  substituted  into  Equation  (100)  gives 


A  -7  * 

u(z)-2a  e  n  s inh  7  ( z  - z ) . 
n  B  n  B  n  F 


(101) 


(102) 


2. 3. 3. 2  Radiating  modes  (k  <k  ) .  In  this  case,  the  eigenfunction  becomes 


n  B 

A 

u  (z)  -  a  cos(rj  z)  +  b  sin(r?  z) 

nB  nB  n  nB  n 

2  2  2 

where  q  -k  -k  .  With  the  boundary  condition  at  z-z  we  get 

n  B  n  F 

b  -  -  a  cot (n  z  ) 
nB  nB  n  F 

which  substituted  into  Equation  (103)  yields 

A  3. 


,  .  nB 

u  (z)  -  — — ; - r 

nB  Sin  (n  Z  ) 

n  F 


sin  rj  (z  -z)  . 
n  F 


(103) 


(104) 


(105) 


Now  we  may  rewrite  Equations  (102)  and  (105)  as 


f 


u  (z)  ■  a 

nB  nB 


sinh  7  (z-z),  k  >k  (trapped) 
n  F  n  B 


(106) 


sin  f?^(zF-z),  k<kfi  (radiating) 

Now  we  use  the  boundary  conditions  at  the  top  of  the  basement  to  match 
Equations  (93)  with  Equations  (106).  From  the  continuity  of  pressure  we  get 


p  u  (z  )  -  p  a 

J  nJ  B  B  nB 


sinh  7  (z  - z  ) ,  k  >k 
n  F  B  n  B 

sin  n  (z  -z  ) ,  k  <k 

n  F  B  n  B 


(107) 


and  from  the  continuity  of  the  normal  particle  velocity  we  obtain 


a  u'  (z  )  —  a  < 

3  nJ  B  nB 1 


7  cosh  7  (z  -z  ) ,  k>k 

n  n  F  B  n  B 

r)  cos  t)  (z  -z  ) ,  k<k 

n  n  F  B  n  B 


(108) 


and  dividing  the  two  equations  we  get  the  characteristic  equation 

7  coth(7  D), 


a  u'  (z  )  1 

V(k)  •  An-  E-  +  ~  s 

PJ  U  ( Z  )  PB 


n  J  B 


r)  cot  (t]  D)  , 

n  n 


k  >k 

n  B 


k  <k 

n  B 


(109) 


where  D  ■  z  -  z  . 

F  B 
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To  avoid  singularities  in  searching  for  the  eigenfunctions  and 
eigenvalues  lets  rewrite  Equation  (109)  for  the  trapped  modes  as 


A  A 

W  (k)  -  —  u  (z  )  +  —  u’  (z  )  tanh(7  D) 
T  ttnJB  pnJB  n 

J  J 


and  for  the  radiating  modes  as 


W_(k)  "  7T  u  ,(2J  cos(»?  D)  +  — -  u'  (z  )  sin(»?  D) 

R  CL  nJ  B  n  p  n  J  B  n 

J  J 


(110) 


(111) 


which  is  still  the  function  of  interest  since  we  are  only  interested  in 

finding  the  zeroes  of  the  characteristic  equation.  If  the  trail  wavenumber 

equals  an  eigenvalue,  the  characteristic  equation  for  that  mode  becomes  null. 

To  search  for  these  eigenvalues,  the  locally  convergent  Newton's  method  is 

6  6 

used  to  converge  to  the  zero  of  the  characteristic  equation. 


From  Equation  (107)  we  get 


a  -  —  u  (z  ) 

nB  p  nJ  B 
B 


csch  7  (z  -  z  )  ,  k  >k 

n  F  B  n  B 


(112) 


CSC  rj  (z  -  z  )  , 
n  F  B 


k  <k 

n  B 


therefore  we  obtain  for  the  eigenfunction  at  the  basement  the  expression 


c  sinh  7  (z  -z) 

_ n  F 

a  a  sinh(7  D)  ’ 

u  (z)  -  p  /p  u  (z  )•<  .  ,  . 

nB  J  B  nJ  B  Sin  ff  (z  -z) 

_ n  F 

sin  rj  D  ’ 

^  n 


k  >k 

n  B 


k  <  k 

n  B 


(113) 


To  eliminate  any  upward-reflected  waves  that  could  have  been  taken  care 
of  by  a  complex  wave -number  we  rewrite  the  eigenfunctions  as 


exp  [7  (z  -z)J 

_ n  F  U  >w 

2  sinh (7  D)  ’  n  b 

n 


U  (z)  -  P  / p  U  (z  M  . .  .  . 

nB  J  B  nJ  B  exp  [  tr)  (z^-z) 

2tsin  r?  D 

v  n 


(lift) 


k  <  k 

n  B 


where  the  up-going  waves  have  been  explicitly  erased  from  the  equations.  If 


we  define 
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N  pt  * 

r  n  J 

F  ~  —t -  u  (z  ) 

n  2 p  nj  B 
B 

the  normalized  eigenfunction  of  the  trapped  mode  becomes 

r 


U  (z)  -  — — r— - — 

nB  sinh(7  D) 


exp [7  (z  -z)] 

n  r 


and  the  radiating  mode  has  a  real  and  an  imaginary  part  given  by 

r 

Re(u  )  -  7— 7 — ;rr  sin[y  (z  -  z)l 

nB  sm  (r?  D)  n  F  J 


(115) 


(116) 


(117) 


and 


r 


Im(u  ) 
nB 


n 

sin(»?  D) 


cos [7  (z  - 

n  F 


Z)  ] . 


(118) 


n 

Since  this  model  of  an  oceanic  waveguide  contains  no  absorption 
whatsoever,  it  is  easy  to  visualize,  from  Equations  (113)  and  (117),  that  the 
radiating  mode  contains  half  of  its  energy  propagating  downwards  and  the 
other  half  is  propagating  upwards.  The  imaginary  part  in  Equation  (118) 
represents  the  absorptive  effect  from  the  basement  when  we  artificially 
erased  the  up-going  wave  from  the  solution.  However,  absorption  is  present 
in  all  layers  of  the  media  and  we  will  include  it  in  Section  2. A. 


2.3.4  Up-Layer  Matching  Algorithm 

When  the  sound  speed  is  strongly  upward-refracting  over  many  of  the 

deeper  waveguide  layers,  it  was  found  that  better  numerical  stability  was 

afforded  by  stopping  the  downward-propagated  solution  at  some  layer  D  and 

propagate  the  solution  upward  from  the  basement  to  layer  U-D+l  with  the  final 

Wronskian  relation  coded  to  match  solutions  at  the  D/U  layer  interface  zy 

rather  than  at  z  . 

B 

At  the  resilient  "false"  bottom  the  eigenfunction  is  automatically  set 

A 

-25 

to  the  small  value  u  (z  )  -  10  for  the  radiating  modes  since  it  is  an 

nB  B 
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exponential  function.  For  propagating  the  solution  upwards,  the  boundary 


conditions  at  zg  give: 


p  u  (z  )  -  p  u  (z  ) 
J  nJ  B  B  nB  B 


(119a) 


U'  (z  )  -  7  u  (z  )  (119b) 

J  nJ  B  n  nB  B 

where  f  is  undefined  in  the  isospeed  basement  and  it  is  possible  to  write 

“u<z»>  -  kjv  <120> 

and  after  the  substitution  of  Equation  (73),  the  boundary  conditions  become 

A 

a  Ai  +  b  Bi  -  R  u  (z  )  (121a) 

nJ  JB  nJ  JB  UJ  nB  B 

and 


where 


a  Ai'  +  b  Bi'  -  T  u'  (z  ) 

nJ  JB  nJ  JB  UJ  nB  B 


R  ■  p  /p  , 
uj  j+r  j  J 

T  ■  a  /a,  j-D . J-l 

Uj  j+l  j  J 


(121b) 

(122a) 

(122b) 


T  -  7  /a  . 

UJ  n  J 


(122c) 


Solving  the  two  equations  for  a  and  b  ,  as  done  before,  in  matrix 

nJ  nJ 


form  yields 


R  Bi'  u  (z  )  -  T  Bi  u'  (z  ) 

UJ  JB  nB  8  UJ  JB  nB  B 

A  A 

T  Ai  u'(z)-R  Ai'  u  (z) 

UJ  JB  nB  B  UJ  JB  nB  B 


(123) 


and  generalizing  this  solution  we  have 


a 

R 

Bi' 

A 

u 

(z 

) 

-  T 

Bi 

A 

u’ 

(2 

nj 

Uj 

j  j*i 

D  1 

j+  1 

u  j 

j  j  +  i 

nj+1 

J+  1 

b 

—  n 

T 

Ai 

A 

u' 

(z 

) 

-  R 

Ai’ 

A 

u 

(z 

Uj 

j  j  + 1 

nj+1 

j+ 1 

Uj 

j  J  ♦  i 

n  j  ♦  1 

j+  1 

uj  uj  JJ-i  iij-j  jui  Uj  J  J  1  iij-l  j-J  /10/\ 

^  A  A  (124) 

b  T  Ai  u'  (z  )  -  R  Ai'  u  (z  ) 

njj  Uj  j  j  +  1  nj+1  j+l  Uj  jj  +  1  nj+1  j+l 

which  is  the  recurrence  relation  for  upward  calculations.  In  the  case  of  a 
very  large  positive  value  of  the  argument  of  the  Airy  function  the  matrix  to 


solve  is 
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u  (Z) 

n  J 

A 

-  n 

u'  (Z) 
nj 

11 


12 


E  E 
21  22 


R  u  (z  ) 

Uj  nj+1  j+1 

A 

T  u'  (z  ) 

Uj  nj+1  j+1 


where 


E  (f  )  -  Bi:  Ai(f  )  -  Ai'  Bi(f  .) 

11  nj  J  J  +  1  nj  j  j  +  1  n J 


n  1 


E,  (f 

)  -  Ai 

Bi (f  )  - 

Bi 

12  n 

j  j  j  + 1 

n  j 

j 

j  +  i 

,  (f 

)  -  Bi'  Ai ' (f  )  - 

Ai' 

2  1  n  j 

j  J  +  i 

nj 

j 

j  +  i 

,,(C 

)  -  Ai  Bi' (f  )  - 

Bi 

22  nj 

J  J  +  i 

n  j 

$ 

J  +  i 

-  and 

down- layer 

solutions 

at 

2  . 
U 

n  j 


(125) 

(126a) 

(126b) 

(126c) 

(126d) 


boundary  conditions  are  applied  at  j-U: 


u  (z  )  -  R  u  (z  ) 

nD  U  UD  nU  U 

A  A 

u'  (z  )  -  T  u  (z  ) 

nD  U  UD  nU  U 


which  are  combined  to  obtain 


u  (z  )  u'  (z  )  - 

nD  U  nD  U 


R  u  (z)u'(z) 

UD  nU  U  nD  U 

A  A 

T  u'(z)u  (z) 

UD  nU  U  nD  U 


or  equating  them  to  the  equation 


R  u(z)u'(z)-T  u(z)u'(z) 

UD  nU  U  nD  U  UD  nD  U  nU  U 


(127a) 

(127b) 

(128) 

(129) 


which  is  satisfied  only  if  the  trail  value  is  an  eigenvalue  of  the  acoustic 

waveguide.  Therefore,  the  up- layer  Wronskian  or  characteristic  equation  is 

defined  as 

A  A  A  A 

®(k)  -  R  u  (z  )  u'  (z  )  -  T  u  (z  )  u'  (z  )  (130) 

UD  nU  U  nD  U  UD  nD  U  nU  U  ' 

which  is  zero  under  the  condition  just  mentioned. 
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2.4  ABSORPTION  AS  A  PERTURBATION 

Now  we  will  consider  the  complex  eigenequation, 

,2 


dz 


-  a  (z)  +  [*2(z)  -  k2}  u  (z)  -  0 

2  n  n  n 


(131) 


where  we  redefine  the  wave-number  as  A(z)  -  k(z)  +  iea(z) ,  t  is  used  here  to 
keep  track  of  the  effects  of  every  term  in  the  resulting  approximate  complex 
eigenequation  and  it  will  be  set  to  unity  at  the  end  of  the  calculations, 
a(z)  is  the  attenuation  coefficient  in  nepers/meter,  and  k(z)  -  w/c(z).  The 
complex  wavenumber  in  Equation  (131)  makes  the  eigenvalues  and  eigenfunctions 
complex.  If  q(z)  «  k(z),  then  we  can  use  the  perturbation  method  to  obtain 
a  more  accurate  transmission  loss.  In  this  case  we  will  write 


and 


\  „  CO)  (1)  ^  2  (2) 

a  (z)  - »  u  +  e  u  +  £  u 

n  n  n  n 


x  CO)  ,  .(1)  2.(2) 

A  +  £  A  +  £  A 

n  n  n 


(132) 


(133) 


which  substituted  in  the  complex  eigenequation  gives, 


r  .2 

d  N.o;  w  X  ,  ,  2  2..  .(0)  ,(1)  2.(2) 

-  +k  (z)+2t£k(z)a(z)  -  £  a  (z)-A  -£A  -e  A 

I  .  2  n  n  n 


(0)  (1)  2  (2) 
U  +  £  U  +  £  U 

n  n  n 


(134) 

which  is  an  approximation  to  the  complex  eigenequation  Equation  (131)  due  to 
the  expansions  Equations  (132),  (133). 

Combining  the  £°  terms  of  this  equation  gives  the  0th  order  solution  to 


the  problem,  or 


d  (0)  .,2.  .  ,<o>.  (0)  A 

-  u  +  k  (z)  -  A  ]  u  -  0 

,  2  n  n  n 

dz 


(135) 


which  is  the  unperturbed  eigenequation  that  has  been  solved  for  the  purely 


real  eigenvalues  A 


(0)  «  k2  and  eigenfunctions  u<0) 
n  n  n 


u  .  This  unperturbed 
n 
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eigenequation  corresponds  to  Equation  (45). 

Combining  the  terms  with  t 1 ,  which  corresponds  to  the  first  order 

perturbation  terms,  gives 

,2 

d  (1)  ,  .  .<1>,  (0)  ,  ,,2,  .  .(0),  (1)  . 

-  u  +  [2ck(z)a(z)  -  A  ]  u  +  [k  (z)  -  A  ]  u  -0  (136) 

.  2  n  n  n  n  n 

dz 

where  the  unperturbed  eigenfunctions  are  normalized  by  Equation  (46) ,  which 
in  the  new  notation  becomes 


J 


p(z)  u  (z)  u  (z)  dz  -6 

n  in  nm 


(137) 


where  z^  is  the  depth  of  the  resilient  bottom  of  the  basement. 


Multiplying  Equation  (136)  by  pu<0>  and  integrating  yields 
-bdu  _  b  .  b 

f  <0)  n  ,  ,  f  <0)ro.,  .  .  .  .  ,(1),  (0)  f  (0)r,2.  ,<0),  (1) 

p u  -  dz  +  pu  [2'-r:(z)Q(z) -A  Ju  dz  +  pu  [k  (z)-A  ]u  dz-0 

Jn,2  Jn  nn  Jn  nn 


dz 


(138) 


where  using  -he  orthonormality  condition  of  the  unperturbed  eigenfunctions  in 
the  second  term  of  this  equation,  integrating  by  parts  twice  the  first  term, 
and  using  the  boundary  conditions  at  every  interface  to  cancel  out  the 

surface  contributions  gives 

*  .2  (°)  z  z 

-bdu  .  b  _b 

f  (1)  n  .  0.  f  (0),  /  \  /  \  .  f  ( 0 )  r .  2  (0).  (1),  (1) 

pu  -  dz  +  2t  pu  k(z)a(z)u  dz  +  pu  [k  (z)-A  Ju  dz  -  A 

Jn,2  Jn  n  Jn  nn  n 

0  dz  0  0 


(139) 

and  with  the  help  of  Equation  (135)  the  first  and  second  integrals  cancel  out 
giving  us  the  expression 

Z 

b 

A(1)-  21  |  p  k(z)  q(z)  I  u(0)  1 2dz  (140) 

n  J  n 

0 

which  is  the  first  order  perturbation  term  for  the  eigenvalue  and  its  values 
are  purely  imaginary. 

Now  we  write  the  perturbed  part  of  the  eigenfunction  under  the  basis  of 
the  unperturbed  part  since  this  is  an  orthonormal  basis,  i.e. 
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(i)  r 

u  "  / 
n  Li 

m 

and  substitute  in  Equation  (136)  to  obtain 

I*, 


A  u 

nm  m 


(0) 


(  141) 


,2  (  0  ) 
d  u 

m  ,  2  .  -<0).  (0) 

— —  +  (k  (z) -  A  )  u 
2  n  n 


dz 

(0) 


+  [2tk(z)a(z)-A(1)]u<0)  -  0  (142) 

n  n 


then  multiply  by  pu^  and  integrate  as  done  before.  Integration  by  parts 
twice  cancels  a  few  terms,  and  the  orthonormality  condition  yields 


a  2c  f  ,  .  .  ,  .  c°)  co)  , 

A  -  -  pk(z)a(z)u  u  dz 

nl  .  (  0  >  ,  <  0  )  J  '  n  1 


(143) 


A  -  A  o 
1 


which  is  in  terms  of  the  unperturbed  eigenfunctions  and  eigenvalues,  is 
directly  proportional  to  the  absorption  coefficient,  and  it  is  a  purely 
imaginary  term. 

In  the  cases  of  trapped  modes,  where  the  imaginary  part  of  the 
eigenvalues  is  extremely  small,  we  can  rely  on  the  rapid  convergence  of  the 
perturbation  method  and  forget  about  a  second  order  perturbation  term.  When 
radiating  modes  are  taken  into  account,  we  must  make  sure  that  the  imaginary 
part  of  the  eigenvalue  is  much  smaller  than  the  real  part.  We  may  stop  the 
eigenvalue  search  when  the  real  part  of  the  eigenvalue  is  about  50  times 
larger  than  the  imaginary  part,  since  the  higher  modes  will  not  propagate 
very  far.  Hamilton  has  estimated  that  the  shear  attenuation  coefficient  is 
about  200  times  higher  than  the  compressional  attenuation  coefficient,  and 
the  shear  speed  is  from  2  to  30  times  smaller  than  the  compressional 

42  44  47 

speed.  '  ’  However,  Werby  and  Tango  found  that  the  shear  to 

4  9 

compressional  attenuation  coefficient  is  only  of  the  order  of  two.  Also, 


Beebe  and  Holland  found  this  ratio  to  be  of  the  order  of  six. 


3  4 


Even  though 


there  is  some  disagreement  about  the  attenuation  coefficients,  it  is 
suspected,  from  the  ratios  given  here,  that  the  imaginary  part  of  the  shear 
wavenumber  may  not  be  as  small  as  the  compressional  one  is.  In  this  case,  it 
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is  feasible  to  attempt  to  solve  the  complex  eigenequation  without  the  first 

order  perturbation  method  or  try  to  get  a  second  order  perturbation  term  and 

find  out  how  large  these  effects  really  are. 

The  approach  here  is  to  continue  with  the  second  order  perturbation 

which  is  needed  mostly  by  the  trapped  modes  when  shear  waves  are  included  and 

2 

by  the  radiating  modes.  The  e  terms  of  Equation  (136)  into  a  second  order 

equation  gives 

,2 

d  (2)  .,2.  .<0).  (2)  ....  .  .  .  .  ,(0).  (1)  .  2.  .  (2) .  (0) 

— -  u  +  [k  (z)-A  ] u  +  [2tk(z)a(z)-A  lu  -  [a  (z)+A  ]u 

2  n  nn  nn  nti 

dz 


(144) 


which  multiplied  by  pu(0)  and  integrated  as  done  with  the  first  order 

n 

eigenvalue  leads  us  to  the  equation 

i  b  *  1 

.(2)  - .  f  .  .  (0)  (1),  .<l)f  (0)  (1).  f  2.  .  |  (0 )  i  2 

A  -  2t  pk(z)o(z)u  u  dz-A  pu  u  dz  -  pa  (z)  u  dz 
n  J  nn  njnn  J  1  n  ' 

0  0  0 

(145) 

where  substituting  Equation  (141)  and  the  orthonormality  condition  of  the 
unperturbed  eigenfunctions  gives 

z  z 

A<2)-  2tV  A  rpk(z)a(z)u(0,ut0)dz  -  fpo2(z)  |  u<0)  |  Zdz  (146) 

n  Li  m  J  nm  J  'n 

m  0  0 

or  with  Equation  (143)  we  get  the  simpler  form 


.  (2)  V  .2  ,.<0)  .(0).  fb  2,  .  ,  (0 )  |  2  , 

A  -  )  A  (A  -  A  )  -  pet  (z)  u  dz 
n  L*  nm  r  m  J  n  ' 


(147) 


which  is  purely  real  and  a  much  smaller  term  since  it  is  proportional  to  a2. 
If  we  write 

(2)  T  _  (0>  .  0. 

u  -  )  B  u  (148) 

n  Lt  nm  m 
m 

then  Equation  (144)  becomes 
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E* 

m 


nm  ,  2 

dz 


(ojr  „  M  2  ,<o)  (0)  r*  .  uk  to)  .  i  ,(2),  (o> 

u  +)  B  [k  -A  ]u  +)  A  [2tka-A  ]u  -  [a  -A  ]u 

m  nm  n  m  L»  nm  nm  nn 


(149) 


which  multiplied  by  pu^0)  and  integrated  using  integration  by  parts  and  the 
orthonormality  condition  reduces  the  equation  to 

(A<0)-  A<0>)  B  -  A(1)A  -  V  A  A  , l*n  (150) 

1  n  nlnnlLnmlm 

m 

2 

which  makes  B  purely  real  and  directly  proportional  to  a  . 

We  have  already  assumed  layers  of  constant  density  in  order  to  simplify 
the  elastic  wave  equation.  Therefore,  we  may  define  an  element  of  a  G-matrix 
as 


G  - 

nm 


i  -  1 


P. 

J 


.j+  1 

k  (z)  a  (z)  uC0)(z)  ut0>(z)  dz  -  G 

J  j  j  nj  mj  mn 

z 

j 


and  that  of  an  H-vector  as 

H  ■ 


j  + 1 


z 

j+l 


j*i  z  j 


(151a) 


(151b) 


Note  that  all  elements  of  the  G-matrix  are  purely  imaginary  and  symmetric, 
while  those  of  the  H-vector  are  purely  real.  These  integrals  must  be 
evaluated  in  order  to  calculate  the  perturbed  parts  of  the  eigenvalues  and 
eigenfunctions . 

Now  the  first  order  perturbation  term  of  the  eigenvalue,  Equation  (140), 
becomes 


A(1>-  G  (152) 

n  nn 

which  tells  us  that  the  diagonal  components  of  the  G-matrix  are  the 
first-order  perturbation  term  of  the  eigenvalues.  The  second  order  term, 
Equation  (145),  simplifies  to 
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r - -  G 

Au>-  y  — — —  -  h 

n  (_ _ _  (  0  )  (  0  )  n 


(153) 


A 

l^n  n 


which  substituted  into 

-2  .(0)  Cl)  .(2)  /1C/N 

*  “  X  +  X  +  X  (154) 

n  n  n  n 

gives  the  perturbed  eigenvalues  of  the  problem.  Since  A(1)  is  the  only 

n 

contributor  to  the  imaginary  part  of  the  eigenvalue,  we  may  define 

X2  -  A<0,+  A(2>  (155) 

n  n  n 

as  the  real  part  of  the  eigenvalue.  Then  to  obtain  k  from  Equation  (154)  we 


expand  its  square  root  as  follows: 

c  l ) 


c  l ) 


k  =  X 
n  n 


r  A'"  , 

1/2 

r  A  -> 

1  +  n 

2  X 

tl 

1  +  n 

X  T 

X2  ' 

2  X2  - 

-  X  + 


(  1  ) 


2  X 


(156) 


2  ( 1 ) 

where  we  have  assumed  that  X  »  X 

n  n 

eigenvalue  will  be  defined  as 


Now  the  imaginary  part  of  the 


( l ) 


2  X 


(157) 


which  is  the  same  expression  in  page  20  of  Reference  24  or  in  Equation  (474) 
of  Reference  14  where  only  first  order  perturbation  has  been  used.  By  the 


same  token,  the  real  part  is  given  from  Equation  (155) 

a(2> 

X~/A<0>+  " 


2/1 


( o  ) 

n 


(158) 


where  it  is  assumed  that  At2>  «  A<0>  and  the  same  power  expansion  has  been 

n  n 

used.  Equation  (157)  is  a  crude  approximations  made  by  many,  however,  and 
they  are  avoided  by  taking  the  complex  square  root  of  Equation  (154). 

As  the  first  order  correction  of  the  eigenfunction,  Equation  (143) 
simply  becomes 
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n  1 


Anr  ;(o)_  Aio, 

n  1 


,  n^l 


and  for  the  second  order  correction,  Equation  (150),  we  get 

-  G  G  r — »  G  G 

_ nn  n  1  +  \  _ n  m  ml _ 

nl“  (AC0,-'a(0>)2  ^  (\<0)-  A<0,)(A<0,-  A(0)) 

n  1  m^n  n  m  n  1 

which  substituted  into  the  equation 


a  (z)  -  u(0)+  y  [A  +  B  ]  u(0> 

n  n  L _ ,  nl  n  1  1 


n^l 


(159) 


(160) 


(161) 


gives  a  better  estimate  of  the  eigenfunction.  These  perturbation  terms  are 
expected  to  shift  the  eigenvalues  into  the  first  quadrant  of  the  complex  plan 
as  displayed  in  Figure  6  where  the  maximum  value  of  the  real  part  of  the 
eigenvalues  is  calculated  using  the  minimum  speed  in  the  sound  speed  profile. 

It  is  left  to  properly  evaluate  of  the  G-matrix  and  the  H-vector  in 
Equations  (151).  Instead  of  numerically  integrating  these  equations,  it  is 
possible  to  make  use  of  some  properties  of  the  functions  contained  in  the 
integrands . 

In  every  layer  we  may  write  c  (z)-ct  +(z-z^  )g  where  the  subscript  T 
stands  for  the  value  at  the  top  of  the  layer.  Now  the  real  part  of  the 
wavenumber  may  be  expanded  assuming  speed  gradients  smaller  than  unity  as 
follows , 


u> 


k  (z)-w/c  (z)-  - 
J  3  c  -g  .  Z 

Tj  Tj 

u> 


1  + 


c  -  g  z 

Tj  6J  Tj 


c  -g  z 
Tj  6j  Tj 


1  - 


c  -  g  z 

Tj  &J  Tj 


(162) 


and  if  we  define  c  ■  c  -g  z  ,  K  ■  w/c  ,  and  M  ■  K  g  /c  ,  nen  we  have 

0j  Tj  bj  Tj’  Tj  Oj  j  Tj°j'  Oj 


k  (z) 
j 


k  -  n  z 

Tj  j 


(163) 

which  is  linear  with  depth  if  the  sound  speed  gradient  is  smaller  than  unity. 
The  attenuation  coefficient  in  Equations  (151)  is  greatly  dependent  on 
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63 


the  frequency  of  the  source  as  shown  by  Urick  in  the  equation 


0.1  f  40  f 

a  -  -  +  - 


+  2 . 75xl0"4  f2 


(164) 


1  +  f  4100  +  f 

where  this  absorption  coefficient  is  in  decibels  per  kiloyards  and  the 
frequency  is  in  kilohertz.  The  dependence  upon  depth  is  given  by  the 
expression 

a  -  a  (1-1 . 93xl03z) 
o 


(165) 


where  a  is  the  attenuation  coefficient  at  the  surface  of  the  ocean.  This 
o 

depth  dependence  is  so  small  that  we  will  take  it  to  be  constant  along  every 
layer.  In  this  case  we  may  take  the  attenuation  coefficient  outside  the 
integrals  and  leave  it  inside  the  summation  across  the  layers. 

The  eigenfunctions  at  all  the  layers,  with  the  exception  of  the 
basement,  are  given  by  a  linear  combination  of  Airy  functions  as  given  in 
Equation  (73).  To  perform  the  integrations  in  Equations  (151a  and  b)  with 
this  property  of  the  Airy  functions  we  will  use  Gordon's  formulas  given  as 

_  57,67-68 

follows , 

A[a(z+b) ]  B[a(z+b) ]  dz  -  (z+b)AB  -  A'B'/a  (166a) 

JzA[a(z+b)  ]B[a(z+b)  ]dz-AB(z2-zb-2bZ)/3  +  (A'B+AB'  )/6a2  +  (2b-z)A'B'/3a 

(166b) 

A'B  -  AB' 

— r -  (166c) 


j*  A[a(z+bj)]  B[a(z+b2)]  dz 
JzA [ a ( z+b  ^ ) ] B [ a ( z+b  ^ ) ]  dz 


a  (b  -  b  ) 
1  2 


2A’B’ 


a‘(b  -b  )2 
1  2 


a3(b  -b  )2 
1  2 


A'B-AB' 


a  (b  -b  ) 
1  2 


where 


and 


A[a(z+bi) J-a^Ai [a(z+bj) j+b^Bi [a(z+bi) 


AB(bi+bz+2z) 

a3(b  -b  )2 
1  2 

(166d) 


( 166e ) 
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B[a(z+bz) ]-a2Ai[a(z+b2) ]+b2Bi [ a(z+b2) ] 


(166f ) 


represent  any  linear  combination  of  Airy  functions  as  our  eigenfunctions  are. 
Now  we  have  modified  Equation  (151a)  to  the  form 


•J  L 

r - '  pj  ♦  1 

-2c)  pa  [  (K  -M  z)  u(0)[a  (z+b  )  ]  u<0>  [a  (z+b  ) 
n  L ,  J  j  ^  Ij  j  n  j  j  n  j  mj  j  mj 


dz  +  G 


(167) 


where  G  is  the  contribution  from  the  isospeed  basement  to  be  calculated 

nm 

later  on  in  this  section.  This  equation  may  be  rewritten  as 


J 

G  -  2i)  pa  (  K  3 
nm  L _ ,  J  J  TJ 


-  M  3  )  + 

nmj  j  nmj  nm 


(168a) 


where 


•j  + 1 

f  (0)r„ 

u  ia. 

i  nj  J 


(z+b  )]  u  [a  (z+b  )  ]  dz 
nj  mj  j  mj 


(168b) 


„j+ 1 

(0)  r„ 
z  u  [a 

4  "j  J 


(z  +b  )]  u  [a  (z+b  )]dz 
nj  mj  j  mj 


(168c) 


are  the  integrals  to  be  solved  analytically. 

The  diagonal  elements  of  the  3 -matrix  are  solved  substituting  Equation 


(166a)  into  Equation  (168b)  giving 


(z+b  )u2 
n  j  n  j 


(z)  -  u'  (z)/a  J 
n  J  J  Z 


(169a) 


where  the  superscript  (0)  has  been  dropped  for  simplicity.  Equation  (166b) 


is  used  to  evaluate  the  n-m  elements  of  the  5-matrix  giving 

if  "I z 

5  -4  (z2-zb  -2b2  )u2  (z)  +u  (z)u'  (z)/a2  +  (2b  -z)u' 2(z)/a  : 

nnj  J  nj  nj  nj  nj  nj  j  nj  nj  j  Z 


(169b) 

The  off-diagonal  elements  of  the  3-matrix  is  calculated  using  Equation  (166c) 


and  is  simplified  to  the  form 
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3  - 

nmj 


T  2 


nm  j 


j+l 


(169c) 


nm  j  J 


where  JE  ■  a  (t  -  t>  )  and  £  ■  u'  (z)  u  (z)  -  u  (z)  u'  (z).  Finally, 

nmj  j  nj  mj  nmj  nj  mj  nj  mj 

Equation  (166d)  is  used  to  calculate  the  off-diagonal  elements  of  the 
3 -matrix  to  give 


3  -  *'2. 

nmj  nmj 


2u'  u'  +  (zJE  +2 a  JE  1  )£  -  a  (t  +b  +2z)u  u 

nj  mj  nmj  j  nmj  nmj  j  nj  mj  nj  mj 


T  2 


j+l 


(169d) 

where  IE  and  £  become  null  for  n-m  and  Equations  (169a)  and  (169b)  must  be 


used. 


The  normalized  basement  eigenfunction  portion  is  represented  by 

s  inh  7  ( z  -  z ) 

- .  ,  ■ , — =rr — ,  k  >k  (trapped  modes) 

sinh(7  D)  n  B  ' 


(0).  .  .  (0).  .  . 
u  (z)  -  p  /p  u  (z  ){ 
nB  J  B  nJ  B  ‘ 


sin  n  (z  -z) 
n  F 


(170) 


- : — — ,  k  <  k  (radiating  modes) 

I  sin(f?  D)  n  b  &  7 

n 

2  2  2  2  2  2 

where  D-z  -z  is  the  thickness  of  the  basement,  7  -k  -k  ,  and  n  -k  -k  . 

F  B  n  n  B  n  B  n 

For  the  trapped  modes  we  substitute  in  the  equation 


-F 

„(E)  0,  W  f  (0).  .  (0).  .  , 

G  -2tp  a  ——  U  (z)  u  (z)  dz 

nm  B  B  C  J  nB  mB 


(171) 


B  z 


the  first  expression  in  Equation  (170)  to  yield 


,2  r  <  0  >  /  \  ,  2 

tp  a  w  [u  (2  ) J 

G<B1-  — - — - - -  [7  :coth(7  D)  -  D  csch2(7  D)  ] 

nn  PC  n  n 


B  B 


for  n-m  and  it  becomes 


0,2  (  0  )  .  ( 0  >  ,  . 
2 tp  a  u>  u  (z  )  u  (z  ) 
,(B)  J  B  nJ  B  mJ  B 


G  - 


^bcb^K> 


[7  coth(7  D)  -  7  coth(7  D) 
tn  m  n  n 


(172a) 


(172b) 


when  n*m.  Both  equations  simplify  to 

2 tp  q  w  u  (z  )  u  (z  ) 
£<B)^  JB  nJ  B  mJ  B 

nm  ^=c0(t  +T  ) 

E  B  n  m 

when  7  D,  and  7  D  are  much  greater  than  unity. 


(173) 


For  the  radiating  modes  we  must  separate  the  solution  to  the  cases  where 
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n-m  and  where  n*m.  In  the  first  case,  Equations  (170)  and  (171)  yield 


.  2 

tp  a  u 

„(B>  J  B  f  (  0  )  ,  s  !  2 
G  -  -  [u  (z  )  ] 

nn  PC  n  J  B 

B  B 


2  -  1 

esc  (r)  D)  -  r\  cot (17  D)  ,  n-m  (174a) 

n  n  n 


which  corresponds  to  Equation  (172a)  for  the  trapped  modes,  and  for  the 
second  case  we  get 


-  .  2  <  0  )  .  \  <°>/  \ 
2 to  a  u>  u  (z  )  u  (z  ) 
(B)  J  B  nJ  B  mJ  B 


G  - 

nro 


Vb  <V\> 


[ »7  co c(r)  D)  -  r)  cot (r?  D)  ]  (174b) 

mm  n  n 


which  is  compatible  to  Equation  (172b)  for  the  trapped  modes. 

But  this  is  not  all.  We  will  also  have  trapped-radiating  mode 
combination.  In  the  case  where  mode  n  is  trapped  and  mode  m  is  radiating  our 
equations  yield 


rt  :  2  (  0  )  ,  .  <  0  )  . 

2tp  a  <0  u  (z  )  u  (z  ) 
(B)  J  B  nJ  B  mJ  B 


G  - 

nm 


Vb<\+\> 


[7  COth(y  D)  -  T)  cot (r)  D)  ] 
n  n  mm 


(175) 


which  is  like  a  combination  of  Equation  (172b)  and  Equation  (174b).  Since 
the  G-matrix  is  symmetric,  we  do  not  have  to  evaluate  the  integral  where  n  is 
a  radiating  mode  and  m  is  a  trapped  one . 


Finally,  we  also  have  to  evaluate  H  ,  given  by 

n 


•J  M. 

H  -  )  p  a2  I  uZ  fa  (z+b  )]  dz  +  H(B> 

n  [_ jjJ  nj'-j  nj  J  n 


(176) 


which  by  the  use  of  Equation  (166a)  we  get 


H  -  \  p  a2  3  +  H(B> 

n  (_ _ _  j  j  nnj  n 


(177) 


where  3  is  given  in  Equation  (169a)  and  it  is  easy  to  find  that 

nnj 

z 

-F 


„<B)  2  2  .  t  ,  \  ^(B) 

H  «  p  a  u  (z)  dz  -  (c  /w)q  G 

n  B  B  .  nB  B  B  nn 


(178) 
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2.5  THE  NORMALIZATION  FACTOR 


The  unperturbed  eigenfunction  must  satisfy  the  orthonormality  condition 
in  Equation  (46).  They  automatically  satisfy  the  orthogonality  condition 
when  searching  for  the  zeroes  of  the  characteristic  equation  because  this  is 
an  implicit  property  of  the  eigenequation,  Equation  (45).  However,  the 
eigenfunctions  are  not  automatically  normalized  because  a  function  which  is 
any  constant  multiplied  by  Equation  (73)  also  satisfies  Equation  (45). 

Hence,  we  must  solve  Equation  (82)  for  the  only  constant  that  will  normalize 
the  eigenfunctions  in  Equation  (83).  This  orthonormal  eigenfunction  will  be 
the  unperturbed  function  that  satisfies  Equation  (135).  Since  the 
eigenfunction  at  the  isospeed  basement  is  different  from  that  of  the  other 
layers,  we  will  divide  Equation  (82)  into  the  basement  and  "Airy"  layers 
contribution 


N~Z  -  IB  +  X 


(179) 


where 


*  -  P.  f  |a  Ai(f  )  +  b  Bi(f  )  1 2dz 
n  L _ ,  J  J  nj  nj  nj  nj 


(180) 


j  “  i  j 

is  the  contribution  of  all  the  "Airy"  layers,  or  layers  other  than  the 

basement,  and  with  the  help  of  Gordon’s  formula  Equation  (166a)  we  obtain 

j 


E  'j  . 


(z+6  )uZ  +  a  1  u'Z 
nj  nj  j  nj 


J+l 


(181) 


j-  i 


where  a  and  *  are  defined  in  Equation  (84).  The  coefficient  IB  is  the 

J  nj  n 

basement  contribution 


IB 

n 
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and  after  substituting  Equation  (170)  for  the  eigenfunction  and  integrating 
we  get 


IB  - 

n 


2  2  . 
",  u„j(V 
20. 


7  1coth(7  D)  -  Dcsch2(7  D) ,  k  >k 
n  n  n  n  B 

D  csc2(fj  D)  -  q  1  cot (t)  D)  ,  k  <k 

n  n  n  n  B 


(183) 


which  gives  us  the  complete  solution  for  the  normalization  coefficient  in 
Equation  (179). 

A  question  that  may  have  come  to  the  reader's  mind  is  that  of  the 
orthonormalization  of  the  perturbed  eigenfunction.  We  have  ensured  the 
orthonormalization  of  the  unperturbed  eigenfunction.  However,  when  using  the 
perturbation  method,  do  we  automatically  have  orthonormal  perturbed 
eigenfunctions?  Direct  substitution  of  Equation  (161)  into  Equation  (46) 
gives 

f  ★  V  ’  ★  * 

\  p  u  udz-B  +  B  +  )  (A  A  +  B  B  +  A  B  +  A  B  ),  n*<m 

J  n  m  nm  an  ^ _ (  In  lm  In  lm  In  lm  la  In 

l**n  .  tn 


(184) 


and  for  the  "renormalization  coefficient"  we  get 


n 


1  +  L  (i\j2  +  i^i2’ 

m^n 


(185) 


which  proves  that  the  perturbed  eigenfunctions  are  not  orthonormal  unless  the 

perturbation  coefficients  in  Equations  (159)  and  (160)  are  zero.  However,  we 

2 

have  shown  that  A  and  B  are  of  the  order  of  a  and  a  ,  respectively,  which 
means  that  the  perturbed  eigenfunctions  are  very  close  to  orthonormal.  The 
reason  for  this  behavior  of  the  perturbed  eigenfunctions  is  that  the 
perturbation  method  is  just  an  approximation  to  the  problem  of  solving  the 
complex  eigenequation  with  complex  solutions  and  boundary  conditions. 

However,  we  have  found  that  this  correction  to  the  normalization  factor 
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affects  the  seventh  significant  digit  of  the  perturbed  eigenfunctions  and  we 
may  simply  assume  the  perturbed  eigenfunctions  to  be  orthonormal.  The 
renormalization  factor,  Equation  (185),  will  be  calculated  for  this  report 
though,  but  there  is  no  simple  way  to  make  the  perturbed  eigenfunctions 
purely  orthogonal.  The  Schmidt  orthogonalization  procedure69  is  extremely 
time  consuming  for  the  small  correction  to  be  made.  Equation  (184)  will 
simply  be  used  to  find  out  how  much  error  is  introduced  in  the  calculation 
when  the  perturbation  method  is  used  instead  of  solving  the  complex 
eigenequation  directly. 
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2.6  MODE -COUPLING  FOR  RANGE  DEPENDENCE 


Now  that  the  depth  part  of  the  solution  has  been  found,  it  is  next  to 
treat  the  range  part  of  the  solution.  Consider  the  inhomogeneous  Helmholtz 
equation. 


[V2  +  k2(r  ,z)  ]  <p(r  ,z)  - 


2mr  5(r>  5<Z-V 


(186) 


where  the  solution  will  be  written  in  the  almost-separated  form 

H 


<p( r.z)  -  )  *  (r)  u  (r,z) 

L _ ,  n  n 


(187) 


n  *  1 


which  embodies  the  "adiabatic  range -variation"  method.  The  unknown  is  4>  (r) 

n 

and  the  function  u  (r,z)  is  the  nth  range- independent  eigenfunction  at  the 

n 

range  segment  where  r  is  contained.  If  u  (r,z)  are  taken  as  the  "basis" 

Ti 

eigenfunctions  that  satisfy  the  equation 


— -  u  (r.z)  +  [k2(r , z)  -  k2(r)]  u  (r,z)  -  0 
.2  n  n  n 

3z 


(188) 


where 


J 


p(z)  u  (r,z)  u  (r,z)  dz  -  6 

n  m  run 


(189) 


then  the  inhomogeneous  Helmholtz  equation  becomes 
N.  (  . 

[(V2  +  k2)*  lu  +  2  V  f  •  V  u  +  *  v2ul--~  S  (r)  S(  z-z  ) 

Lr  nnJn  rn  rn  nrn  Z7 XV  0 

I 

n*  1 


n 

L 


where 


3r 


i  d_ 

2  +  r  3r  ' 


(190) 


(191) 


Multiply  both  sides  by  p( z)  u  (0,z)  and  integrate  to  obtain  the  range 

JT5  “ 

equation 
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d  ^  1  d  x  ,2,  , 

— r  +  -  +  k  (r) 


L  dr 


r  dr 


where 


and 


dtf 


1  *  _ 

*  <r>  -  -  ^  Hr)  p(z  )  u  (0,z  )  -  2)  -r^  W  (r) 

n  <djrr  o  n  o  u  ar  nm 

m 

-  Tf  (r)M"  (r) 
Lt  i d  nm 


M'  (r)  ■  f F  p(z)  u* ( r , z )  f-  u  (r.z)  dz 

nm  J  m  OX  n 


fZF  * 

[a2  l  a 

J  p(z)  um(r,z)  | 

L  ar2  +  ?  ar  J 

u  ( r , z )  dz 

n 


(192) 


(193a) 


(193b) 


or 


where 


V  (r) 

nm 


M"  (r)  (r)  +  V  (r) 

nm  r  nm  nm 


zf  *  d 2 

p(z)  u  (r.z)  — -  u  (r , z)  dz 

m  3r2  n 


(193c) 


(193d) 


are  the  coupling  coefficients  of  the  inhomogeneous  equation  and  these  are 
null  when  the  eigenfunctions  have  no  range  dependence.  Note  that  Equation 
(192)  is  a  set  of  N  equations  for  (r)  coupled  in  all  the  depth 

n 

eigenfunctions.  This  is  the  equation  to  be  solved  by  the  coupled  normal-mode 
method.  To  simplify  this  equation  it  is  desirable  to  eliminate  the  first 
derivative  by  introducing  a  new  function  defined  by 


#  (r)  -  r'1/Z  f  (r) 

n  n 


( 194) 


which  transforms  Equation  (192)  to  the  form 

d2  1  2  . 

— r  +  — ~  +  k  (r) 

i  .  2  »  2  n 

Ldr  4r 


*  df 

f  (r) - 4—  6( r)  p( z  )  u  (O.z  )  -  2V  M'  (r)  ~  - 

n  0  1  /  2  On  0  L  nm  ar 


2n  r 


Y  V  (r)  f  (r) 

L  nm  m 


(195) 


or  all  N  coupled  equations  can  be  written  in  the  matrix  form 
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”  +  K2(r) 
.  2  .2 
dr  4r 


F(r)  -  U(r)  -  2  M' (r)  F' (r)  -  V(r)  F(r)  (196) 


where 


F(r)  - 


f:(r) 

f2(r) 


f  (r) 

N 


(197) 


and  its  derivative  is  the  derivative  of  every  element  in  the  matrix.  Also  we 


»ve 


K2(r) 


k  (r)  0 

l 

0  k2(r) 

2 


and 


U(r) 


-1 


2wr 


6 (r)  p(z  ) 
1/2  r  o 


0 

0 


• •  k  (r) 

N 


u  (0 , z  ) 

1  0 

u, (°>zn) 

2  0 


u  (0,  z  ) 

N  0  J 


(198) 


(199) 


which  helps  in  the  simplification  of  the  range  equation. 

The  method  consists  in  dividing  the  range  -  dependent  environment  into 
range -independent  segments  where  the  first  segment  is  the  only  one  with  a 
source  and  the  last  one  is  semi-infinite  as  qualitatively  shown  in  Figure  7. 
The  coupling  will  occur  in  the  other  segments  and  the  radial  boundary 
conditions  are  easily  satisfied  as  long  as  the  slope  of  the  bottom  is  small. 
Under  these  assumptions  we  can  simplify  the  coupled  range  equation  for  every 
segment  and  find  its  solutions. 
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2.6.1  Range  Segment  Number  1 

In  this  segment  we  assume  range  independence  and  note  that  this  is  the 
only  segment  containing  the  omnidirectional  source.  Under  these  conditions 


Equation  (196)  decouples  to  form 


- —  +  —  +  K2(r  ) 
Ldr2  4r2  0 


F(1)(r)  -  U(r) 


(200) 


which  has  the  homogeneous  solution 


,<1). 

f  (r) 

n 


1/2 


(1)  u(l>  ..  0  (1)  (2)  .  0  . 

a  H  (k  r)  +  B  H  (k  r) 

n  0  n  n  0  n 


(201) 


where  k  -  k  (r  ),  and  the  alphas  and  betas  are  coefficients  to  be  determined 

n  n  0 

via  the  radial  boundary  conditions.  The  inhomogeneous  term  U(r)  will  be 
taken  into  account  later  on  when  matching  the  field  to  the  source.  For 
simplicity,  this  solution  can  be  rewritten  as 


f(1>(r)  -  (r)  2H  (r)]  a 

n  v  r\  ^  J 


(1) 

n 


(202) 


where 


-*(i  1 

a  — 

n 


(i) 

t 

n 

(i) 


and 


JH  (r)  -  r1/2  H*j>(k°r)  ,  j-1,2. 

n  On 


(203) 


(204) 


2.6.2  Range  Segment  Number  M+l 

In  this  segment  we  also  assume  range  independence,  but  there  also  is  no 
source.  This  makes  Equation  (196)  uncoupled  and  homogeneous.  Also,  due  to 
its  semi- inf inite  property,  we  will  only  accept  the  radiation  solution  to  the 

equation 
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^  +  -L_  +  K2(r  )' 
,  2  .  2  M 

dr  4r 


Ftmi>(r)  -  0 


which  is 


-(M+l)  .  .  1/2  (M+l)  „(1>  M  . 

f  (r)  -  r  ct  H  (k  r) 

n  n  0  r. 


(205) 


(206) 


where  k  -  k  (r  ).  For  simplicity  we  can  write  these  functions  as  elements 

n  n  M 

of  the  matrix  equation 


where 


F(M+1)(r)  -  A  H(r) 


..  ,  .  1/2  „(1)  ..  M  . 

H  (r)  ■  r  H  (k  r). 

n  On 


(207) 


(208) 


2.6.3  Range  Segments  Number  2  <  i  <  M 

These  segments  are  range  dependent  as  a  whole  and  the  coupling  terms  are 
kept.  However,  the  source  term  is  zero  and  to  simplify  the  equation  we 
assume  these  segments  to  be  at  a  range  much  larger  than  a  wavelength.  Under 
these  assumptions  the  range  equation  for  these  segments  become 


r~~2  +  k2<V 

Ldr2 


Fa,(r)  +  2  M'  (r  )  F'U>(r)  +  V(r  )  F(i>(r)  -  0  (209) 

i  i 


where  r  -  (r  +  r  )/2  and  the  matrices  are  constant  values  averaged  to  be 

i  i-1  i 

determined  at  the  boundary  between  every  range  segment.  A  series  of 
transformations  is  next  performed  in  order  to  repeatedly  diagonalize  matrices 
and  regroup  terms  until  the  modes  are  uncoupled. 

To  solve  Equation  (209)  we  define  as  the  range  -  independent  matrix 

that  diagonalizes  M'(r  )  to  have  diagonal  elements  X  as  follows, 

i  n 


A  =  s'V  S 
1  1 


r  A  ^  ...  O', 
0.  .  .  A 


(210) 
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where  S  S  -1  is  a  useful  property  satisfied  by  this  transformation  matrix. 
If  we  multiply  the  left  hand  side  of  the  range  equation  by  S^1  then  this 


—  +  X2  +  V  I  r<r)  +  2  A  T' (r)  -  0 


equation  becomes 


where  X2  -  S'VCr^S^  V  -  S^V  A  -  S'Vs^  and  T(r)  -  sVi>(r), 
Now  we  will  define 

a  (r)  -  exp(-A  r) 

n  n 


(211) 


(212) 


as  the  elements  of  the  diagonal  matrix  o( r) ,  where  a  a  —  0  is  satisfied,  and 
also  define 


r(r)  -  ff(r)  7(r) 


(213) 


which  substituted  in  the  equation  gives 


r  g2  2  1 

I  ~  +  X2  +  V  a  7  +  2  A  (a  7'  +  o' 7)  -  0 


(214) 


and  now  multiply  by  a  to  the  left  to  obtain 


- —  +  T  7(r)  -  0 
*2 


(215) 


where 


T  .  a'\  XZ+  V  )  a  -  A2 
-  a'1  S'1(K2+V)S  o  -  A2 

l  l 


(216) 


7(r)  -  o'V1  F(i)(r) 


(217) 


have  simplified  the  equation  considerably,  but  they  are  still  coupled  because 
T  is  not  diagonal.  They  are  finally  uncoupled  by  finding  the  matrix  S  that 


diagonalizes  T  to  have  elements  A  with  the  transformation 


A  -  S  T  S 
2  2 


A  .  .  .  0 
.  1 

0.  .  .A 


(218) 


where  S  S'1-!!  is  satisfied.  Multiplying  the  left  hand  side  of  Equation  (215) 
2  2 
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by  takes  us  to  the  final  form 


+  A 


L  dr 


f(r)  -  0 


which  has  the  solutions 


(219) 


?(r)  -  s'S(r) 


q^cos (q^r)  +  ^isin(qir) 
o^cosCq^r)  +  ^sinCq^r) 


(220) 


a  cos(q  r)  +  /3  sin(q  r) 

N  N  n  nN 

1/2 

where  q  -  A  ,  and  a  ,  /9  are  the  unknown  to  be  determined  by  the 
n  n  n  n 

application  of  the  boundary  conditions  and  all  superscript  (i)  indicating  the 
range  segment  of  evaluation  has  been  momentarily  dropped  for  clarity.  Taking 
all  the  transformations  that  connect  F(r)  to  f(r)  it  is  easy  to  find  that 


F(r)  -  [S i<r(r)  SJ  ?(r)  -  V  P(r) 
or  back  to  the  original  form  we  get 


n 


(r) 


N 


m  *  1 


(i  )  /  (i  )  N  ,  ^(i  )  .  ,  (i  )  N  ' 

a  cos(q  r)  +  sin(q  r) 

m  ra  m  m 


and 


(221) 


(222) 


FU)(r)  -  U1  7,i( r) 


(223) 


where  n  -  1,2,...,N  is  the  mode  index,  i  -  2,3, ...,M  is  the  segment  index,  T1 
is  not  a  diagonal  matrix,  and  the  unknown  coefficients  are  in  P1 (r) .  Now 
that  the  radial  functions  have  been  found,  it  is  next  to  match  them  to 
satisfy  the  radial  boundary  conditions . 
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2.7  THE  RADIAL  BOUNDARY  CONDITIONS 

The  radial  boundary  conditions  are  the  same  in  principle  to  those  for 
the  depth  functions.  The  only  difference  is  in  the  coordinate  used  to 
satisfy  for  the  continuity  of  pressure  and  normal  particle  velocity.  Since 
the  range  segments  are  small  and  a  linear  fit  is  used  to  connect  the  nth 
eigenfunction  radially,  then  it  is  assumed  that  the  eigenfunction  and  its 
radial  derivative  are  radially  continuous.  Substituting  Equations  (194), 
(187),  and  (28)  into  Equation  (79)  gives 

f<U(r  )  -  f U+1) (r  )  (224) 

n  i  n  i 

for  the  continuity  of  pressure  in  the  radial  direction  at  the  interface  r  . 

The  next  boundary  condition  is  given  by  taking  the  radial  derivative  of 
the  velocity  potential 


d(p 

dx 


v— -  r  d#  du  ^ 

-  )  |— —  u  +  # 

/  .  I  dr  n  n  dx 


(225) 


n  *  1 


where  the  last  three  functions  in  the  parenthesis  are  already  continuous  and 


the  first  function  is  the  one  to  be  satisfied  by  writing 


df 

r 

d F 


(  i  ) 


df 

T 

dF 


( i+i ) 


(226) 


where  i  -  1,2, . . . ,M. 


Apply  the  boundary  conditions  at  r  to  get 

M 

A  H(r  )  -  ?M(r  ) 

M  M 


and 


where 


A  H'  (r  )  -  Ti”  (r  ) 

M  M 


,M  -1/2  (1)  ,  1  -3/2„(l  M 

kr  H  ’(kr)-^r  H  (k  r  ) 

nM  0  nM  2  0  nM 


(227a) 


(227b) 


(228) 
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and 


ajM,  .  ,  M  M  .  .  M  .  ,  aH  M  .  M  . 

(rM)  "  'Qn  Sln(q„rM)  +  **„ 

nM  n  n  nn  nn  nn 


(229) 


To  further  simplify  the  calculations,  let's  define  a  hyper- space  vector  as 

r 


«u ) 

a  — 


♦(i) 

t 

1 

*U) 

t 

2 

*(1) 

i 

H 


(230) 


and 


Cj(r  ) 


which  has  the  inverse  form 


( 0  (r  )  ] 

n  i 


and  the  hyper -space  matrix 

Cj(r  ) 


cos(qJr  )  sin(qJr  ) 


■qJsin(qJr  ) 

n  n  1 


cos(qjr  ) 

n  i 

sin(qJr  ) 


qj  cos (qj  r  ) 


(231) 


-sin(qJr  )/qJ 

n  i  n 

cos (qJr  )/qJ 


n  i  n 


(232) 


C3(r  ) 
1  1 


0 


C2(ri> 


(233) 


with  an  inverse  given  by  the  inverse  of  the  individual  diagonal  elements  in 
this  matrix.  Finally  it  is  helpful  to  define 


f  H  (r  ) 

I  H' (r  ) 

n  M 


,  n-1,2 . N 


(234) 


and  the  hyper-space  vector 


fl 


(235) 
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then  the  matching  of  the  radial  functions  at  r  becomes , 


and 


N 

<M)  ,  (M)  .  ,  a(M)  .  .  (M)  .  Y  ’  I"  -1  1  „  .  . 

a  cos (q  r  )  +  &  sin(q  r  )  -  )  [U  J  A  H  (r  )  (236) 

n  nM  n  n  M  /  ^  y  mra 

- -  u  J  nm 

m*  1 


n 

(M)  (M)  .  ,  (M)  x  ,  ^(M)  (M)  .  (M)  x  \  1  [  r77CMK-lA  1  ,  N 

•a'  q  sm(q  r  )  +  0  q  cos(q  r  )  -  )  U  A  H' (r  ) 

„  „  n  M  »  »  n  M  [  J  ™  “  “ 


ro  ■  1 


or  in  hyper-space  matrix  form,  both  equations  combine  into  the  equation 


<r«„>  =r  -  E  [  Cu'H'r'A ; 


H 


(237) 


(238) 


0-1 


or 


(239) 


S(M)«  [C<M)(r  )]-1[U(M>]-1  A  fi 

M 

which  is  a  recurrence  relation  between  the  last  semi -  inf inite  segment  and  the 
last  finite  segment  coefficients  to  be  determined. 

For  the  intermediate  boundaries,  i-2 , 3 , . . . ,M- 1 ,  the  functions  and  their 


derivatives  are  matched  to  obtain 


and 


?u)(r  )  -  [U<i)]’1  U(i+1)  fti+1)(r  ) 

i  i 


Tu)'  (r  )  -  [UCi) ] _1  U(i+1)  (r  ) 

i  i 


(240) 


(241) 


th 


Combining  the  n  row  of  each  and  gathering  both  sets  of  equations  leads  to 

N 

C(i)(r  )  a(i)  -  'T  [U(i) ]  1  U(i+1)C(i  +  1,(r  )  (242) 

n  in  _ i  m  1  m 


m  ■  1 

and  shifting  to  hyper-space  gives 

a  -  C  (r)]  [U  0  C  (r  )  a  ,  i-2 , 3 , . . . ,M- 1 

i  i 

(243) 

which  is  the  recurrence  relation  for  the  unknown  coefficients  at  the  finite 
range  segments. 
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After  applying  the  boundary  conditions  of  the  remaining  interface  and 

writing  both  sets  into  one  matrix  equation  we  get 

K 


K  a(1)-  r  U(2,C(2)(r  )  o(2) 

n  n  L _ ,  nm  m  1  m 


(244) 


where 


0*1 


H  (r  ) 

n  1 

XH'  (r  ) 

n  1 


H  (r  ) 

n  1 


H'  (r  ) 

n  1  J 


and  in  hyper- space  notation  it  is  found  the  third  recurrence  relation 


where 


»(1>  „-l  (2)  (2),  «*(2> 

a  -  K  u  C  ( r  )  a 


K'1  0 

l 


(245) 


(246) 


(247) 


N  J 

is  the  matrix  that  uncoupled  the  source  segment  coefficients,  q(1)  and  /?(1), 

n  n 

( 2  )  ( 2  )  (M+ 1 ) 

thus  providing  them  in  terms  of  a  and  /9  ,  or  ultimately,  a  .  After 

n  n  n 

multiple  substitution  of  the  calculated  hyper- space  equation  we  get  a 
relationship  between  the  unknown  coefficients  in  the  source  segment  and  the 
unknown  coefficient  in  the  semi- infinite  segment  given  by, 


a1-K*1U2C2(r  )  fl  |  [  C1  ( r . )  ] _1  [U1  ]  "V+1Ci+1(r  )i[CM(r  )  ]"1(UM]‘1A  fl 

1  i-21  1  11  ” 


(248) 

where  the  parenthesis  at  the  superscript  have  been  momentarily  eliminated  to 
write  the  equation  in  one  line.  This  relationship  can  be  greatly  simplified 
to  the  form 

S{1)-  X  A  fl  (249) 

where 
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X  -  K'1  f]  |  U{i>Cli>(r  )  [C(i>(r  )]‘X  [Uu>]*1|. 

l-z*-  1  1 


th 


The  n  hyper-row  is  given  by 


1 


X  A  H 

TUB  SB  D 


(250) 


(251) 


where  A  -  a  ,  and  back  out  fro.n  the  matrix  notation  the  coefficients  for 

ffl  m 

the  first  range  segment  are  given  by 

Q(1)-  y  A  (  Xu  H  +  X12  H'  ) 

n  La  m  run  m  no  m 


and 


$U)-  y  A  (  X21  H  +  X22  H'  ) 

n  La  m  nro  so  run  id 


(252) 


(253) 


where 

■  xu  xi2  ■ 

X  -  ""  nm  (254) 

nm  x2i  x22 

L  nro  nm  J 

and  rM  has  been  dropped. 

It  is  now  time  to  take  into  consideration  the  behavior  of  the  source  at 
r-0.  This  is  the  last  matching  to  be  performed  to  obtain  the  necessary 
number  of  equations  to  evaluate  all  the  unknowns  in  a  closed  form.  In  this 
first  segment  it  was  found  that  the  homogeneous  solution  to  the  inhomogeneous 
and  uncoupled  equation  is 

(255) 


,  ,  .  (  1  )„(  1  )  0  _  fl<  1  )„<2  )  0  . 

(r)  -  q  H  (k  r)  +  0  H  (k  r) 

n  n  0  n  n  0  n 


and  in  the  limit  as  k  r  — >  0  we  find  that 

n 


which  gives 


1  )  /,  0  ,.(2  )  ..  0  2i  .  ,,  0  . 

H  (k  r)  -  -  H  (k  r)  — >  —  ln(k  r) 

On  On  7T  n 


.  21  ( l  >  .( l ) .  ,  /,0  s 

't  (r)  — )  —  (q  -  f)  )  ln(k  r)  . 
n  7T  n  n  n 


(256) 


(257) 


An  asymptotic  solution  can  also  be  calculated  by  integrating  the  uncoupled, 
inhomogeneous  equation  over  a  small  cylinder  containing  the  source,  and 
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taking  the  limit  as  the  radius  goes  to  zero.  The  procedure  gives 


d2* 


„  „  .  «.  d*  - 

j  n  dr  +  |  i  -r-ndr  +  |  k2(r)  #  (r)  dr  - 

J  j_2  J  r  dr  J  n  n 


Pi O  u  (0 ,  z  ) 


0  n 


o  dr 


where  the  second  term  is  integrated  by  parts  to  yield 


2* 


f  -JSll  dr 

oJ  r 

(258) 


d*  ,*  *  ,* 


dr 


+  —  +  f  f  r'2+  k2(r)]  *  (r)  dr  -  - 

r  J  ^  n  J  n 


P ( z  )  U  (0 ,  z  ) 


0  n 


*  6(r) 


1  o  o 


2n 


dr 


(259) 


and  in  the  limit  as  a  — )  0  only  the  slope  at  r-0  and  the  integral  over  the 


delta  function  remains  providing  the  expression 

d*  (r)  p(z  )  u*(0,z  ) 

n  0  n  0 


dr 


2nr 


(260) 


or 

*  (r)  ->  -  tX-  p(z  )  u*(0 , z  )  ln(k°r)  (261) 

n  i. IT  0  n  0  n 

which  can  now  be  equated  to  the  other  asymptotic  solution  providing 

Q(1).  0(1)-  l/h  p(z  )  u*(0 , z  )  -  V  (262) 

n  n  0  n  0  n 

which  is  the  relationship  between  q(1>  and  /3C1).  Substitution  of  the 

n  n 

calculated  coefficients  into  this  equation  gives 

V  -  V  [(X11  -  X21)  H  +  (X12  -  X22)  H'  ]  A  ■  y  B  A  (263) 

n  nm  run  m  nm  run  tn  cn  run  m 

m  m 

or  in  matrix  form 

V  -  B  A  (264) 

and  after  inverting  the  matrix  B  (if  it  is  not  singular)  the  first  unknowns 
are  obtained  with  the  equation 

A  -  B'JV.  (265) 

After  this  matrix  is  calculated,  it  follows  to  calculate  the  unknowns  of 
the  source  segment  with  Equations  (252)  and  (262)  .  The  coefficients  of  the 
other  segments  are  calculated  using  Equations  (239)  and  (243).  Finally,  the 
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range  functions  are  obtained  with  Equations  (201),  (206),  and  (222)  where  the 
Hankel  functions  are  calculated  with  the  complex  eigenvalues  obtained  in 
Section  2.U. 

Note  that  this  coupling  method  is  in  a  closed  form  that  requires  no 
iteration  that  may  slow-down  the  computation  time.  It  is  left  to  evaluate 
the  coupling  coefficients  which  is  done  in  the  next  section. 
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2.8  REDUCING  THE  COUPLING  INTEGRALS 


The  coupling  integrals  are  given  in  Equations  (193)  where  the 
eigenfunctions  satisfy  Equations  (188)  and  (189).  Under  the  assumptions  of 
an  environment  that  slowly  varies  with  range,  it  is  possible  to  use  Gordon's 
formulas  in  Equations  (166)  to  solve  the  coupling  integrals. 

The  range -dependent  environment  is  divided  into  range  - independent 
segments  as  shown  in  Figure  7.  The  eigenvalues  and  eigenfunctions  have  been 
calculated  for  the  middle  of  the  segments  at  r  .  A  method  of  evaluation  is 

fZj.  *  Q 

to  multiply  Equation  (188)  by  the  operator  dz  p(z)  u  (r,z)  —  to  yield 

oj  m  dr 

the  expression 


rZF  *  a 

J  ^  ai 


a_ 

dr 


dz  Un 


^  fZF  *f3k2'j  rZF  *  2  3un 

dz  +  J  puJarhdz  +  J  k  (r,z)  ^  dz  - 


2 

rzv  *  rak  1  rzv  *  3u 

pu  T^udz+  pu  kZ(r)  J-?-  dz  (266) 
J  m  lor  In  J  m  n  or 


where  n*m.  Since  the  eigenvalue  is  independent  of  depth,  the  first  integral 
in  the  right  hand  side  of  this  equation  vanishes  due  to  the  orthogonality 
condition.  Integrating  by  parts  twice  the  first  integral  in  the  left  hand 
side  of  Equation  (266)  gives 


j  + 1 

D 


*  fL 

r  3u 

n  j 

m j  dr 

dz 

3u  3u  z 

m  n 


dz  dr 


z  a  u 


■ 

J  Z  0 


+  k  (r , z) 


dz 


fZF  * 

i  2> 

3k 

j  p  U  U 

J  m  n 

,  3r 

★ 

u  dz  - 

m 


3u 


(267) 

where  no  approximation  has  been  made  yet.  The  complex  eigenfunctions  and 
eigenvalues  are  used  in  Equation  . 195)  to  find  the  complex  f  (r) .  The 

n 

coupling  matrices  are  supposed  to  be  complex  also,  but  under  the  condition  of 
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a  slowly  varying  range  dependence  the  real  part  of  the  coupling  integrals  is 
very  small  compared  to  the  other  terms  in  Equation  (196)  and  the  imaginary 
part  is  even  smaller  because  the  imaginary  part  of  the  eigenfunctions  and 
eigenvalues  is  much  smaller  than  the  real  part.  Assuming  that  the  imaginary 
part  of  the  eigenfunction  is  much  smaller  than  the  real  part,  Equation  (188) 
may  be  substituted  in  the  integral  at  the  left  hand  side  of  Equation  (267)  to 
get  the  approximate  expression 

M'  (r)  -  (Mv  (r)  +  M*  (r))/(k2(r)  -  k2(r)),  n*m  (268) 

ntn  ran  ran  n  m 


where  the  volume  contribution  is  given  by 

j  +  i 


M  (r) 

ran 


and  the  surface  contribution  is 

J+ 1 


\  r j+: 

)  P.  \  U  u 
/ _ ,  J  J  _  rnj  nj 


r3k2. 


3r 


dz 


(269a) 


M°  (r). 

nm 


j-l 

j 

r  3u 

3u 

du  -| 

b 

tnj  bt 

n  j 

IS 

n 

3z 

3  z 

3r 

j  +  l 


(269b) 


j  - 1  j 

where  all  imaginary  parts  are  neglected  and  the  stratified  medium  of 


isodensity  layers  is  used. 

Since  the  range- independent  wavenumber  squared  has  been  written  in  the 


form 


k2(z)  «  k2(z  )  + 
j  j  j 


k2 (z  ) -k2 ( z  ) 
j  j + 1  j  j 


z  -  z 
j  +  l  j 


(z  -  z  ) 

J 


then,  by  linear  interpolation,  we  obtain 

2 


3k  (r,z) 

j _ 

dr 


=»  ft  z  +  B 
r  ji  ji 
1 


where 


(k2  -  k2  ) - (k2  -  k2  ) 

^  ^  j  -*■  1  ,  1  j  ,  i _ j  +  1  ,  i  -  1  j.l-l 

ji*  (Az)  (Ar) 

j  i 

i  *  (k2  -  k2  )/(Ar )  -  ft  z  , 

Ji  j .  i  j.i'l  i  Ji  J 


(270) 


(271a) 


(271b) 

(271c) 
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Ur). 


r  -  r 

i  i-1 


and 


(Az)  ■»  z  -  z  , 
j  j+i  j 


k2  ■  k2  ( r  ,  z  ) 

n,m  j  m  n 


(271d) 

(271e) 

(271f) 


Now  the  volume  contribution  is  given  by, 

j  +  1 


Mv  (r  )  *  V~V  [  J+1u  (r  , z )  u  (r  z)  [k  z  +  IB  ]  dz  (272) 

nm  i  l_ _ ,  J  J  nj  i  mj  i  ji  ji 


j-1 


or 


Mv  (r  )  -  )  p  (A  3  +  IB  ©  ) 

nm  i  l (  j  j  i  nm  j  i  j  i  nm  j  i 


(273) 


j-i 

where  3  and  ©  are  Airy  integrals  of  the  same  form  as  the  ones  in 
nmj i  nm j i 

Equations  (168)  .  Since  the  imaginary  part  of  the  eigenfunctions  have  been 

neglected,  there  is  no  sense  in  using  the  second  order  perturbation  term  for 

the  evaluation  of  these  integrals.  Also,  the  summation  will  go  up  to  j-J 

because  it  has  been  assumed  a  range -independent  isovelocity  basement.  This 

2 

assumption  makes  3kj+i/8r  -  0  and  it  poses  no  restriction  since  we  can 
stratify  the  bottom  to  the  depth  of  interest,  and  consider  from  there  on  the 
so  called  basement  layer  of  large  thickness.  With  the  help  of  Gordon's 
formulas  we  obtain  3  -  3  from  Equation  (169d),  and  6  -  3  from 

nmji  nmj  nmj  i  nmj 

Equation  (169c)  where  nym  and  the  subscript  i  only  indicates  the  range 
segment  into  consideration. 

The  pressure  and  normal  particle  velocity  boundary  conditions  are  useful 


in  reducing  the  surface  contribution  to  the  form 


M  (r  ) 

nm  1 


r3u  3u 

m  r 

dz  dr 


Z  (r  ) 
F  i 


(274) 


and  under  the  assumption  of  an  environment  that  slowly  varies  with  range  this 
contribution  is  negligible  compared  to  the  volume  contribution.  Rutherford 


7  0 


and  Hawker  proved  that,  in  an  environment  with  range  variable  boundaries 
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and  interfaces,  neglicting  the  surface  contribution  results  in  the  violation 
of  the  conservation  of  energy  in  the  radial  direction.  However,  the  acoustic 
properties  in  the  model  in  Figure  7  are  the  only  range -dependent  parameters 
in  the  problem  and  energy  conservation  is  not  violated. 

In  the  case  where  n-m,  Equation  (193a)  may  be  written  as 


*  i  j,  J  h  u»j<r'z)  dz 


(275) 


where  the  imaginary  part  of  the  eigenfunction  has  been  neglected  as  done  for 


the  n*m  case.  Since  the  orthonormality  condition  gives 


H  dz  -  J  I?” dz  +  f 1 \  h  dz  ■ 0 


(276) 


then  Equation  (275)  becomes 


M'  (r)  «  -  |  y. p,  [u2  (r,z)  J* 

nn  2  l _ ,  2  l  nj  dr  J 


(277) 


which  is  another  negligible  contribution  compared  to  the  off-diagonal 
elements  of  the  coupling  matrix  when  the  environment  has  a  weak  range 
dependence . 

In  conclusion,  this  coupling  matrix  is  given  by 


M'  (r  )  « 

nm  i 


M  (r  ) 

hi  i 

k2(r  )  -  k2(r  ) 

n  i  mi 


(278) 


where  only  the  term  with  the  greatest  contribution  has  been  kept.  Note  from 
Equation  (269a)  and  (278)  that  M'  -  -M'  or  that  the  M' (r)  matrix  is 

nm  on 

antisymmetric.  Since  the  eigenfunctions  are  considered  real,  then  the 
elements  of  M' (r)  are  real  too.  Therefore,  the  matrix  is  antihermitian  and 
it  becomes  Herraitian  when  it  is  divided  by  i.  This  way  it  is  easier  and 
faster  to  diagonalize  the  matrix  and  after  the  diagonalization  of  the 
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Hermitian  matrix,  the  diagonal  elements  must  be  multiplied  by  i  to  obtain  the 
eigenvalues  of  M' (r) . 

It  is  next  to  solve  the  coupling  matrix  V(r)  with  elements  given  by  the 
integral  in  Equation  (193d) .  The  solution  is  obtained  with  the  help  of  the 
expression 

f-  u  (r , z)  -  ^  '  M'  (r)  u  (r,z)  (279) 

p^n 

which  is  proved  by  direct  substitution  in  Equation  (193a)  and  the  use  of  the 
orthonormality  condition  in  Equation  (189).  Substitution  of  this  same 
expression  into  Equation  (193d) ,  taking  the  radial  derivative  of  both  terms 
in  the  summation  of  Equation  (279),  and  substituting  Equations  (189)  and 
(193a)  gives 

V  (r)  -  ~  M'  (r)  +  'S  '  M'  (r)  M'  (r)  ,  m*n  (280) 

nm  of  nm  /_ _ t  np  pro 

p^n  ,  ro 

which  is  completely  in  terms  of  the  previously  calculated  coupling  matrix  in 
Equation  (278).  The  radial  derivative  of  the  coupling  element  can  be 
simplified  by  acknowledging  that  the  eigenvalues  vary  faster  with  range  than 
the  volume  contribution  of  the  coupling  matrix.  Then,  at  every  range  segment 


If  n-m,  the  coupling  term  M'  (r)  vanishes  and  the  radial  derivative  in 

nn 

Equation  (280)  is  null.  Since  M'(r)  is  an  antihermitian  matrix,  Equation 
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(281)  becomes 


n 

V  (r  )  -  -  r  [M'  (r  )] ' 

nn  i  l_ _ t  pn  i  J 


(283) 


p- 1 


where  the  squared  of  M' (r)  must  be  very  small,  but  the  summation  of  these 
numbers  makes  V(r)  of  considerable  contribution. 

The  most  suitable  calculation  of  these  coupling  matrices  has  been  a 
source  of  a  large  amount  of  papers.  An  accurate  calculation  of  the  coupling 
matrices  takes  a  finely  spaced  sampling  of  the  eigenfunctions  and  the  range 
segments,  since  radial  derivatives  and  depth  integrations  are  required.  In 
order  to  take  less  computer  memory  and  time,  it  is  appropriate  to  assume  an 
environment  that  slowly  changes  with  range.  Another  way  to  attack  the 
problem  is  to  calculate  the  coupling  matrix  in  Equation  (193a)  direct!’ 

The  range -dependent  environment  is  divided  into  range  -  independent 
segments  as  shown  in  Figure  7.  The  eigenvalues  and  eigenfunctions  have  been 
calculated  for  the  middle  of  the  segments  at  r  and,  for  simplicity,  the 
range  between  boundaries  is  the  same.  If  the  environment  slowly  varies  with 


range,  then  Equation  (193a)  at  the  range  R  becomes 


M'  (R  )=[*  Yp( z) 

nm  1  J 


z  ru(r  ,z)+u(r,zKru(r  ,z)-u(r,z) 

mi  +  l  mi  ii-.-li 


n  i  +  1 


r  -  r 

i  +  1  i 


dz  (284) 


where  the  definition  (AR)  r  -r  and  the  orthonormality  condition  further 

i  i+1  l  J 

reduces  the  equation  to  the  form 


Mnm ^ Ri ) ”  2 (AR) 


J  +  1  Z 

k  xk  r 


★  * 
u  u  -  u  u 
mi  ni+1  roi+1  ni 


dz 


(283) 


j 1 1  j 

where  the  second  subscript  of  the  complex  eigenfunctions  stands  for  the  range 
segmentation  index  in  the  argument  of  the  eigenfunctions  in  Equation  (284)  . 
This  equation  already  gives  a  property  of  the  M'  -matrix  when  the  complex 
conjugate  is  taken.  The  property  is  that 
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M'  -  -  M'*  (286) 

urn  n  m 

and  only  half  of  the  off-diagonal  elements  have  to  be  calculated.  Setting 
m-n  in  Equation  (286)  also  proves  that  the  diagonal  components  of  the  matrix 
are  purely  imaginary.  Since  the  imaginary  part  of  the  eigenvalues  and 
eigenfunctions  are  much  smaller  than  the  real  part,  the  assumption  of 
Equation  (277)  being  negligible  is  valid.  Also  note  that  this  new  result  for 
M' (r)  is  still  antihermitian  and  can  be  divided  by  i  to  diagonalize  the 
resulting  Hermitian  matrix. 

Compare  the  coupling  matrix  as  calculated  in  Equation  (285)  with  the  one 
in  Equations  (268)  and  (269)  .  Equation  (285)  requires  the  integration  over 
complex  eigenfunctions  at  different  range  segments,  while  the  other  method 
only  needs  to  be  integrated  over  unperturbed  eigenfunctions  at  the  same  range 
segment.  The  approximate  solution  in  Equation  (278)  is  therefore  called  the 
adiabatic  approximation  of  M'  (r) ,  since  the  elimination  of  the  imaginary  part 
of  the  coupling  coefficients  implicitly  assumes  that  no  energy  in  the  form  of 
heat  is  allowed  to  be  transferred  from  one  range  segment  to  another.  The 
integration  in  Equation  (285)  does  take  into  account  the  transfer  of  heat 
from  f ne  segment  to  another  when  integrating  over  modes  of  two  contiguous 
range  segments  and  the  imaginary  part  of  the  coupling  matrix  is  a  clear 
indication  of  the  energy  lost  in  the  ith  segment  due  to  this  unadiabatic 
process.  The  point  is  to  integrate  Equation  (285)  as  analytically  as 
possible,  since  numerical  integration  is  lengthy  and  extremely  unreliable 
when  the  rapid  oscillation  of  the  higher  radiating  modes  are  taken  into 
account . 

The  complex  eigenfunctions  in  Equation  (285)  are  the  solutions  of 
Equation  (131).  For  Airy  layers,  the  solution  is  Equation  (73)  where  the 
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coefficients  and  the  argument  of  the  Airy  functions  are  complex.  The 

argument  is  given  by  Equations  (85)  and  (84)  which  are  now  complex  since  the 

proper  representation  of  the  attenuation  coefficient  is  by  a  complex  sound 

speed.  Since  the  sound  speed  may  vary  from  range  segment  to  another,  the 

coefficient  a  in  Gordon's  formula,  Equation  (166c),  is  different  for  both 

functions  Afa^z+b^)]  and  B[a  (z+b  )].  Therefore,  another  analytical  formula 

is  needed  to  solve  Equation  (285)  when  the  eigenfunctions  are  a  linear 

combination  of  complex  Airy  functions  where  the  complex  version  of  a  in 

j 

Equation  (84c)  is  now  given  by 

V  -C-£o"’  (287) 

and  ft  is  the  same  form  as  in  Equation  (84e)  but  complex.  After  calculating 

n  j 

the  complex  coupling  matrix  M'  (r) ,  V(r)  is  calculated  with  Equations  (281), 
(282),  and  (283). 

If  the  sound  speed  profile  does  not  vary  much  with  range,  an  approximate 
formula  to  use  is 

f  A*[a*(z+b*)  ]  B [ a  (z+b)]  dz  -  JLJL:  (288) 

J  a  a  (bi  -  b2) 

which  must  be  used  twice  to  evaluate  Equation  (285)  . 

Note  that  the  adiabatic  approach  requires  less  computations  and  computer 

memory,  also,  in  most  circumstances,  the  oceanic  acoustic  properties  and 

boundaries  slowly  vary  with  range  and  the  adiabatic  approximation  seems  to  be 

the  most  reasonable  approach.  The  imaginary  part  of  the  eigenvalues  and 

eigenfunctions  are  not  used  to  evaluate  the  coupling  coefficients,  but  they 

are  used  to  obtain  the  range  functions  in  Equations  (201),  (206),  and  (222). 

After  the  range  functions  are  computed,  the  particle  velocity  is  given 

by  Equations  (187)  and  (194).  The  transmission  loss  for  the  receiver  at  any 

depth  and  range  larger  than  six  wavelengths  is  calculated  in  the  next  section. 
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2.9  TRANSMISSION  LOSS  CALCULATION 

When  the  signal- to-noise  ratio  at  a  sonar  system  is  greater  than  a 
calculated  detection  threshold ,  it  is  concluded  that  a  target  is  present.  An 
accurate  calculation  of  this  detection  threshold  is  needed  for  avoiding 
unnecessary  false  alarms.  For  passive  sonar  systems,  where  the  sonar  simply 
listens  to  the  target,  the  detection  threshold  is  given  by  the  formula, 

DT  -  SL  +  DI  -  TL  -  NL  (289) 

where  SL  is  the  projector's  source  level,  DI  is  the  receiver's  directivity 
index,  TL  is  the  media  transmission  loss,  NL  is  the  ambient  noise  level,  and 
DT  is  the  detection  threshold  for  the  sonar. 

For  an  active  sonar,  where  a  sonar  system  or  a  sonobuoy  listens  to 
echoes  bouncing  back  from  the  target,  the  calculation  of  the  noise-limited 
detection  threshold  is  given  by 

DT  -  SL  +  DI  -  TL  -  TL  +  TS  -  NL  (290) 

12 

where  TL^  is  the  transmission  loss  from  the  source  to  the  target,  TL^  is  the 
transmission  loss  from  target  to  the  receiver,  and  TS  is  called  the  target 
strength.  The  noise  level  is  replaced  by  the  reverberation  level  for  the 
reverberation- limited  detection  threshold. 

In  both  cases,  an  accurate  calculation  of  the  transmission  loss  is  of 
paramount  importance.  The  transmission  loss  is  defined  as 

TL  -  -  10  log  I(^’Z)  (291) 

where  1^  is  the  time -averaged  reference  intensity  at  1  meter  from  the  source, 
and  the  magnitude  of  Equation  (33)  provides 

I (r  ,  z)  -  |  p  v  |  (292) 
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where  the  pressure,  p,  is  given  by  Equation  (28)  and  the  scalar  potential  by 
Equation  (62).  Therefore,  we  may  decompose  the  time -dependent  pressure  field 
into  its  modal  contributions , 

p(r , z , t)  -  £  pn(r,z,t)  (293) 

n 

where  the  modal  contribution  for  a  range- independent  environment  is 

p  (r.z.t)  -  p(z  )  p(z)  u*(z  )  u  (z)  h”1  (k  r)  etu>t  (294) 

and  direct  substitution  of  the  pressure  into  the  transmission  loss  equation 
is  known  as  the  "coherent"  loss  since  we  are  taking  into  account  the  phase 
factor  of  every  mode.  The  transmission  loss  requires  the  time  average  of  the 
intensity,  and  for  short  pulses  in  active  sonars  the  ratio  of  energy  flux 
density  from  Equation  (32)  must  be  used  in  Equation  (291).  The  ratio  of 
intensities  can  be  used  if  the  pulse  is  long  enough  to  include  the  most 
significant  multipath  modes.  In  shallow  waters,  the  delay  between  the 
multipath  eigenrays  is  very  small  allowing  a  correct  calculation  of  the 
transmission  loss  using  Equation  (291)  with  shorter  pulses. 

If  the  pulses  are  extremely  short,  then  every  multipath  mode  must  be 
taken  into  account  individually  and  the  simple  effects  of  spherical  spreading 
and  absorption  are  taken  into  account  by  the  approximation 

TL  ~  20  log(R)  +  aR  (295) 

where  R  is  the  distance  travelled  by  the  eigenray  which  may  be  calculated 
using  the  eigenray  method  in  Reference  71  in  meters,  and  a  is  the  attenuation 
coefficient  in  Equation  (164)  where  the  units  must  be  converted  from 
kiloyards  to  meters.  However,  the  constructive  interference  of  the  longer 
pulses  yields  a  higher  signal - to-noise  ratio  making  them  more  useful  in 
active  sonar  systems.  For  these  longer  pulses  and  for  longer  distances  Marsh 
and  Schulkin  have  shown  experimentally  that  the  spreading  becomes 
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cylindrical,  the  number  20  in  Equation  (295)  becomes  10,  and  other  correction 
factors  dependent  on  the  sea  state  and  bottom  properties  must  be  added  to 
Equation  (295) . 

When  the  detailed  interference  effects  are  not  of  interest,  they  may  be 
averaged  out  to  yield  smooth  transmission  loss  curves.  This  is  done  by 
summing  the  individual  modal  energies,  rather  than  their  pressures,  and  the 
result  is  called  "incoherent"  loss.  In  this  case,  we  must  time -average  the 
sum  of  the  squares  of  the  modal  pressures.  For  incoherent  loss,  we  replace 


Equation  (291)  with 


TL  -  -  10 


lo*  32  b ?b  U  lRe(p-)] 


(296) 


where  we  may  take  the  average  of  the  sum  to  be  the  same  as  the  sum  of  the 
average  and  using  the  asymptotic  form  of  the  Hankel  function  in  Equation 
(294)  we  get 


p  (r , z , t)  - 


(8*k  r) 

n 


—  p(z  )  p(z)  U*(z  )  u  (z)  ec(knr-u>t-,r/A)  (297) 

1/2  0  n  0  n 


hence 


([Re<p„)ia>  -  \  lpJ2-  nbr  E  p(zo)  ^(z)  VV  un(zb  (298) 

n 

and  the  incoherent  transmission  loss  for  a  range  -  independent  environment  is 


simply 


TL(r,z)  -  -10  log 


i*pz(  z)  ^ 


U2(z  )  U2(z) 
n  0  n 


(299) 


where  the  phases  of  the  summed  modes  are  immaterial  when  taking  the  magnitude 
squared  of  the  eigenfunctions. 

For  the  range -dependent  environment  we  have  calculated  the  velocity 


potential  in  Equation  (187),  therefore  a  transmission  loss  in  terms  of  the 
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velocity  potential  most  be  found.  The  pressure  is  given  in  terms  of  the 
velocity  potential  in  Equation  (28)  and  the  transmission  loss  in  terms  of  the 
pressure  is  given  by 


TL(r,z) 


-20  log 


P(r,z) 

P  (0,z  ) 

^  o  o  J 


(300) 


where  p^  is  the  acoustic  pressure  at  the  distance  of  one  meter  from  the 
source  and  the  calculation  becomes  incoherent  when  the  pressure  is  replaced 
by  the  rms  pressure.  Using  Equation  (28),  the  rms  pressure  becomes 

p(r ,  z)  -  1//2  w  p (r , z)  <p(  r.z)  (301) 

and  the  same  applies  to  the  rms  pressure  at  one  meter  from  the  source  with 
the  exception  that  the  velocity  potential  may  be  given  by 

V  -  A/(4*r2)  (302) 


since  spherical  spreading  may  be  assumed  near  the  source.  Assuming  a  unitary 


source  strength  and  at  a  distance  of  one  meter  Equation  (302)  becomes  <pq=»1  /4w 


and  the  rms  pressure  at  one  meter  is 

P 


u>  p(0,z  ) 
1  o 

4w 


°  /2 

which  substituted,  with  Equation  (301),  into  Equation  (300)  gives 

f4m  p  (r  ,  z ) 


TL(r,z)  -  -20  log 


p(0,zq) 


<p(r  ,z) 


(303) 


(304) 


where  the  complex  absolute  value  of  the  velocity  potential  in  Equation  (187) 
must  be  taken.  Note  that  the  summation  in  Equation  (299)  is  taken  when  the 
velocity  potential  is  calculated. 

The  effects  f  absorption  are  incorporated  in  Equation  (304)  by  the  use 
of  the  perturbed  eigenvalues,  Equation  (154),  and  eigenfunctions,  Equation 
(161),  when  the  attenuation  coefficient  a(z)  in  nepers/meter  as  a  function  of 
depth  is  given.  Absorption  coefficient  is  mostly  given  in  dB/meter  and  we 
must  convert  these  units.  Outgoing  spherical  pressure  waves  in  an  isospeed 
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media  are  attenuated  as 


p(r,t)  -  pQe 


tw( r/c - t) -or 


(305) 


and  the  intensity  from  Equation  (292)  is  proportional  to  the  squared  of  the 
time-averaged  pressure,  therefore,  the  absorption  contribution  of  the 
transmission  loss  in  Equation  (291)  becomes 

(306) 

and  when  the  reference  is  taken  as  one  meter  we  get 

£  -  [20  log(e) ]  a  (307) 

or 


£(r)  -  -  20  log(e"ar) 


a  -  0.1151  £  (308) 

where  £  is  in  dB/meter. 

Another  way  to  describe  the  absorption  of  the  media  is  by  expressing  the 
loss  in  decibels  per  unit  wavelength.  This  value  is  obtained  by  multiplying 
e  by  the  average  sound  speed  of  the  media  and  dividing  the  result  by  the 
frequency  of  the  source. 

The  perturbation  method  in  Section  2.6  may  also  be  applied  to  the  shear 
wave's  attenuation  coefficient.  The  effects  of  shear  waves  from  the  bottom 
will  be  included  in  the  mode  coupling  method  in  a  latter  report.  Other 
approximations  made  are  the  range  segmentation  of  the  range  -  dependent 
environment  as  shown  in  Figure  7.  But  this  should  not  be  a  great  limitation 
since  it  is  possible  to  increase  the  number  of  segments  until  the  calculated 
transmission  loss  curves  reach  a  state  where  it  stays  unmodified  regardless 
of  the  extra  segmentation.  The  same  occurs  with  the  depth  layering  of  the 
sound  velocity  profile  into  constant  speed  gradients.  The  false  resilient 
bottom  added  to  the  model  for  the  discretization  of  the  continuous  spectrum 
can  vary  in  depth  until  the  transmission  loss  curves  at  the  region  of 
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interest  remain  steady  with  increasing  z^.  The  directivity  of  the  source  and 
receiver  are  not  a  problem  to  be  added  to  the  transmission  loss  calculation 
since  they  are  included  in  the  source  level  and  directivity  index  of  the 
sonar  equations,  respectively.  Also,  the  frequency  spectrum  from  a  source 
may  be  modeled  as  a  discrete  number  of  monochromatic  signals  very  close 
together  and  with  different  amplitudes.  This,  as  well  as  the  shear  waves  and 
three  dimensional  effects,  remain  to  be  included  in  the  model. 
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CHAPTER  3 

RESULTS 

A  simple  environment  will  be  considered  in  order  to  minimize  the 
complexity  of  the  problem  and  to  better  study  the  effects  of  range  dependence 
of  the  ocean  floor  and  the  effects  of  its  attenuation.  This  environment, 

^4-15 

also  studied  by  J. Miller,  is  an  upslope  shallow  ocean  similar  to  the  one 

displayed  in  Figure  7.  The  water  has  a  constant  sound  speed  of  1500.0  m/s 
with  a  density  of  1000.0  kg/m3  and  it  is  200.0  meters  deep  from  the  location 
of  the  source  to  a  range  of  5.0  km.  This  depth  decreases  to  60.0  meters  at  a 
range  of  10.0  km  and  it  remains  flat  from  there  on.  The  bottom  is  one  thick 
sediment  with  a  resilient  boundary  at  1200.0  meters,  a  sound  speed  of  1704.5 

3 

m/s,  a  density  of  1150.0  kg/m  ,  and  an  attenuation  coefficient  of  8.4x10 
nepers/m.  Let's  assume  that  the  source  is  a  six-blade  propeller  from  a  112.0 
meters  deep  submarine  moving  at  a  speed  of  25  knots.  Therefore,  as  a  rule  of 
thumb,  the  propeller  is  generating  a  quasi -monochromatic  continuous  acoustic 
wave  with  a  frequency  of  25.0  hertz.  The  range  dependent  environment  will  be 
divided  into  eight  segments  and  the  bottom's  slope  remains  constant 
throughout  the  range  dependent  region  of  the  environment. 

The  resulting  transmission  loss  as  a  function  of  range  and  depth  is 
given  in  the  three-dimensional  plot  of  Figure  8  where  the  range  is  in 
kilometers  and  the  depth  is  in  meters.  A  total  of  18  modes  have  been  used 
for  this  transmission  loss  calculation.  At  the  position  of  the  source,  three 
trapped  modes  have  been  found.  Note  that  just  below  the  source,  sound  is 
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reflected  back  from  the  resilient  false  bottom  regardless  of  the  high 
attenuation  coefficient  given.  However,  the  effects  of  these  standing  waves 
are  negligible  at  the  range  dependent  region  of  the  environment,  where  the 
interesting  phenomena  occurs.  Only  three  trapped  modes  exist  at  the  range 
independent  segment  where  the  source  is  located.  The  source  has  been  set  to 
112.0  meters  deep  to  be  located  at  the  node  of  the  second  eigenfunction  in 
order  to  excite  only  the  first  and  third  mode.  As  the  bottom  of  the  ocean 
decreases,  the  third  mode  gets  cut-off  from  the  water  and  becomes  a  radiating 
mode  that  propagates  into  the  bottom  sediment.  This  explains  the  low 
transmission  loss  region  under  the  bottom  of  the  range  dependent  segments  of 
this  environment.  The  first  mode  stays  in  the  water  column.  Also  note  that 
there  is  some  energy  in  the  second  mode,  since  the  node  is  not  exactly  at 
112.0  meters  deep.  This  mode  is  cut-off  at  a  range  of  about  9.5  kilometers 
from  the  source.  The  grid  in  this  plot  is  divided  into  squares  of  20.0 
meters  in  depth  and  200.0  meters  in  range. 

If  the  depth  of  the  source  is  set  to  a  depth  of  190.0  meters,  then  a 
second  region  of  low  transmission  loss  is  created,  as  shown  in  Figure  9, 
which  is  caused  by  the  convertion  of  the  second  eigenfunction  to  a  radiating 
mode.  Note  the  refracted  wave  created  bounces  back  from  the  false  bottom, 
but  it  is  damped  rapidly  before  it  reaches  the  water  column.  The  lower 
transmission  loss  of  the  refracted  wave  is  caused  by  the  addition  of  the 
second  mode  and  the  fact  that  the  source  is  only  ten  meters  away  from  the 
ocean  floor. 

If  the  source  were  located  at  the  other  side  of  the  wedge,  i.e.,  at  30.0 
meters  depth  in  the  60.0  meters  deep  semi  -  inf inite  range  segment,  then 
down-slope  sound  propagation  occurs  as  shown  in  Figure  10.  However,  the 
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source  is  exciting  the  only  mode  available  in  this  section  of  the  ocean.  As 
the  sound  propagates  towards  the  deeper  waveguide,  only  the  fundamental  mode 
is  excited  and  the  transmission  loss  plot  displays  no  interesting  feature. 
Since  this  behavior  is  expected,  this  result  is  a  good  indication  that  the 
model  is  performing  correctly. 

So  far,  only  isovelocity  layers  have  been  used.  In  order  to  test  the 
performance  of  the  transmission  loss  calculations  with  layers  of  constant 
gradients  in  the  wavenumber  squared,  we  will  include  in  the  up- slope 
environment  a  sediment  layer  with  z  -  300.0  meters  deep,  as  shown  in 
Figure  7.  This  layer  has  a  constant  sound  speed  gradient  of  0.64  s  1  with 
the  sound  speed  at  z  equal  to  the  one  in  the  basement  layer.  Also  the 
density  and  attenuation  coefficient  of  this  layer  is  equal  to  the  one  of  the 
basement.  The  source  is  set  back  to  112.0  meters  deep  in  order  to  excite  the 
fundamental  and  third  mode  only.  The  resulting  transmission  loss  surface  is 
given  in  Figure  11,  where  only  the  portion  over  800.0  meters  deep  is 
displayed  to  better  observe  the  propagation  phenomena  at  the  sediment  layer. 
Comparing  Figure  11  to  Figure  8,  note  that  part  of  the  energy  from  the  third 
mode  is  propagating  though  the  sediment  and  refracts  back  into  the  water 
column,  hence  the  decrease  in  the  transmission  loss  in  the  60.0  meters  deep 
water  column.  Also,  due  to  the  interaction  with  the  sediment,  some  of  the 
energy  of  the  third  mode  is  being  used  to  excite  the  second  trapped  mode 
which  becomes  a  radiating  mode  as  it  propagates  away  from  the  source.  The 
second  mode  is  not  excited  by  the  source,  but  by  the  range -dependent 
waveguide  itself.  This  phenomena  of  energy  interchange  between  the  modes  can 
only  occur  when  the  model  allows  the  coupling  between  the  modes.  This  is  a 
clear  indication  of  the  high  performance  of  this  coupled  normal  mode  model 
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when  low  frequency  sound  transmission  loss  in  a  range -dependent  environment 
is  needed. 
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CHAPTER  4 
CONCLUSIONS 

The  coupled  normal-mode  model  has  been  developed  with  the  effects  of  the 
attenuation  coefficient  in  all  the  layers  of  the  ocean  as  a  first  and  second 
order  perturbation  to  the  unperturbed  eigenequation.  It  has  been  found  that 
the  first  order  perturbation  term  is  purely  imaginary  and  the  second  order 
perturbation  term  is  purely  real  (which  represents  a  correction  to  the 
unperturbed  eigenvalue  and  eigenfunction).  This  second  order  perturbation 
term  is  absolutely  neccesary  when  the  attenuation  coefficients  are  relatively 
high.  This  is  the  case  when  shear  waves  from  the  bottom  sediments  are 
included.  For  cases  where  the  bottom  elasticity  has  a  very  high  compressional 
and  shear  attenuation  coefficient,  higher  order  perturbation  terms  are 
expected  to  be  neccesary.  In  such  cases,  the  perturbation  approximation  is  no 
longer  feasible  and  a  method  for  searching  the  eigenvalues  in  the  complex 
k-plane  is  neccesary.  The  theory  of  normal  modes  with  the  effects  of  the 
bottom's  elastic  properties  and  a  method  of  searching  for  the  complex 
eigenvalues  is  underway  and  will  be  presented  in  future  publications.  Also, 
the  effects  of  three-dimensional  sound  propagation  using  an  extended  version 
of  the  coupled  normal-mode  model  will  be  modeled. 
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FIGURE  1  ACOUSTIC  RAY  THEORY  IN  THE  DEEP  OCEAN  (ACOUSTIC  NORMAL 
MODES  MUST  BE  USED  IN  THE  SHALLOW  OCEAN) 
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FIGURE  2  DEYRIVATION  MODELS 
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FIGURE  5  THE  GENERAL  STRATIFIED  RANGE  INDEPENDENT  OCEAN- 
ENVIRONMENT  MODEL 
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FIGURE  6.  EXPECTED  MIGRATION  OF  EIGENVALUES  WHEN  COMPRESSIONAL 
ABSORPTION  IS  INTRODUCED 
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FIGURE  7.  THE  RANGE-DEPENDENT  ENVIRONMENT  DIVIDED  INTO  M  +  1  RANGE- 

INDEPENDENT  SEGMENTS  WHERE  THE  SOURCE  IS  AT  THE  FIRST  SEGMENT 
AND  A  VERY  DEEP  RESILIENT  BOTTOM  IS  INCLUDED  TO  DISCRETIZE  THE 
CONTINUOUS  SPECTRUM 
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EXCITE  THE  SECOND  MODE) 
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FIGURE  10  DOWN-SLOPE  TRANSMISSION  LOSS  (SOURCE  DEPTH  IS  30  0  METERS 
AND  ONLY  THE  FUNDAMENTAL  MODE  CAN  BE  EXCITED) 


FIGURE  11  UP  SLOPE  TRANSMISSION  LOSS  (SOURCE  DEPTH  IS  1 12.0  METERS;  A  BOTTOM 
SEDIMENT  WITH  POSITIVE  SOUND  SPEED  GRADIENT  IS  INCLUDED) 
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